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PREFACE 


This  report  examines  the  problem  of  partially  saturated 
moisture  movement  from  a  circular  infiltrometer  to  a  water  table. 
While  actually  three-dimensional,  the  problem  can  be  described  by  a 
two-directional  coordinate  system  by  taking  advantage  of  the  axis  of 
symmetry. 

The  report  begins  with  a  discussion  of  Darcy's  law  and 
functional  equations  which  have  been  proposed  by  a  number  of 
investigators  for  defining  an  effective  hydraulic  conductivity  for 
partially  saturated  soils  for  use  in  Darcy's  law.  Readers  familiar  with 
this  background  may  wish  to  begin  their  reading  with  the  section 
"Inverse  formulation,"  and  only  scan  the  sections  "Influence  of 
Parameters  in  Permeability-Capillary  Pressure  Relationship"  and 
"Relative  Permeability  as  a  Function  of  Saturation." 

The  mathematical  description  of  the  physical  problem 
considers  the  magnitudes  of  the  radial  and  the  axial  coordinates  as 


the  dependent  variables  and  the  potential  function  (hydraulic  head) 
and  Stokes'  stream  function  as  the  independent  variables  or  the 
coordinates  which  define  the  plane  of  the  partial  differential 
equation  boundary  value  problem.  The  resulting  inverse  formulation 
reverse  the  usual  role  of  the  variables,  and  has  several  advantages, 
which  are  made  apparent  in  the  report,  over  the  more  conventional 
formulation. 

A  computer  program  has  been  developed  for  solving  the 
resulting  two  nonlinear  partial  differential  equations  with  their 
associated  boundary  conditions  simultaneously  by  utilizing 
techniques  of  finite  differences.  A  number  of  solutions  are  obtained 
by  varying  parameters  which  describe  the  hydraulic  properties  of  a 
range  of  soil  types.  Analyses  from  the  results  from  these  solutions 
are  presented  in  a  number  of  graphs  toward  the  later  part  of  the 
report  which  permit  one  to  predict  a  number  of  features  of  the 
steady  state  flow  under  varying  conditions. 


ABSTRACT 


Solutions  are  obtained  to  the  problem  of  steady-state  partially 
saturated  infiltration  of  moisture  applied  over  a  horizontal  source 
circle  which  moves  through  homogeneous  soils  toward  a  water  table. 
A  commonly   accepted  relationship  between  relative  permeability 
and  capillary  pressure  has  been  utilized  in  conjunction  with  Darcy's 
law    to    formulate    the   mathematical   model.    The   solutions   have 
utilized  an  inverse  formulation  and  have  been  obtained  by  finite 
difference.  The  inverse  formulation  considers  the  magnitudes  of  the 
cylindrical  coordinates  r  and  z  as  the  dependent  variables  and  the 
potential    function    $     and    Stokes'    stream    function  *\>     as    the 
independent  variables  (i.e.  the  problem  is  solved  for  r  and  z  in  the 

$0  plane).  The  approach  used  for  solving  the  problems  is  practical 
with   modem  digital  computers.  The  computer  output  gives  the  r 


and  z  coordinates  at  each  finite  difference  grid  point.  These  values 
can  readily  be  plotted  in  flownet  form  to  show  the  characteristics  of 
the  flow  pattern  at  a  glance.  From  the  solution  results,  the 
distribution  of  capillary  pressure,  relative  permeability,  or  effective 
saturation  over  any  surface  or  plane  of  interest  can  be  obtained.  The 
solutions  indicate  that  significant  radial  movement  of  moisture  (or 
spreading  effect)  occurs  causing  higher  infiltration  rate  at  the  edge 
of  the  source  circle  than  near  the  center.  The  infiltration  rate  is 
closely  related  to  various  soil  parameters  which  characterize  the 
hydraulic  properties  of  soils.  Also  presented  are  several  distributions 
of  the  relative  permeability  or  effective  saturation  on  the  surface, 
along  the  axis  of  symmetry,  and  on  the  plane  including  the  axis  of 
symmetry  and  how  these  distributions  are  related  to  the  soil 
parameters. 
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INTRODUCTION 


Water  movement  through  soils  and  other  porous  media  has 
been  of  interest  to  mankind  since  early  history.  Approximately  one 
hundred  years  ago,  the  scientific  basis  upon  which  to  predict  water 
movement  in  soils  began  to  develop.  The  most  important  contribu- 
tion to  the  quantitative  study  of  the  groundwater  movement  was 
made  by  Henry  Darcy  who,  in  1856,  observed  the  characteristics  of 
downward  flow  of  water  through  sand  filters  and  postulated  the 
now  famous  basic  linear  law  of  flow  of  water  through  porous  media. 
His  postulate  states  that  the  volume  of  water  crossing  a  unit  area  in 
unit  time  is  directly  proportional  to  the  difference  between  heights 
of  water  in  manometers  terminated  above  and  below  the  sand,  and 
inversely  proportional  to  the  thickness  of  the  sand.  In  other  words, 
the  velocity  is  proportional  to  the  negative  hydraulic  gradient.  This 
relation  has  since  become  known  as  Darcy 's  law.  Darcy's  law  applies 
to  porous  media  that  are  saturated  and  contain  only  a  single 
parameter  that  characterizes  the  media  known  as  the  coefficient  of 
permeability. 

Although  Darcy 's  original  experiment  was  made  on  water 
flowing  vertically  downward  through  a  sand,  Darcy's  law  has  been 
shown  to  be  invariant  with  respect  to  the  direction  of  the  flow  in 
the  earth's  gravity  field.  The  generalization  of  Darcy's  equation  for 
saturated  flow  in  three-dimensional  space  can  be  written  as  v  =  -K 
grad  h  (Muskat,  1937,  and  Hubbert,  1956)  in  whichTis  the  velocity 
vector,  and  h  (a  scalar  quantity)  is  the  total  head  defined  as  the 
height  above  a  standard  elevation  datum  of  the  water  column  in  a 
manometer,  and  K  is  a  coefficient  with  units  of  velocity  depending 
on  the  permeability  of  the  porous  media.  While  Darcy's  law  is  stated 
macroscopically,  the  volume  elements  to  which  the  velocity  and 
pressure  refer  are  supposed  to  contain  a  large  number  of  pores  and 
the  dynamical  variables  are  averages  over  a  large  number  of  pores. 
On  a  small  scale  there  are  large  variations  in  size  of  individual  pores 
in  a  disordered  porous  media  such  as  soil.  The  fluid  motions  within 
the  individual  pores  obey  the  Navier-Stokes  equations  of  motion. 
Darcy's  law  is  really  a  statistical  result  giving  the  empirical 
equivalent  of  the  Navier-Stokes  equations  averaged  over  a  large 
number  of  individual  pores.  Thus,  Darcy's  law  is  a  statistical 
macroscopic  equivalent  of  the  Navier-Stokes  equations  of  motion 
for  the  flow  of  water  through  porous  media  (Hall,  1956,  and 
Hubbert,  1956). 

Darcy's  law  itself  does  not  permit  the  solution  of  a  particular 
problem  of  flow  of  water  through  porous  media.  Rather  it  is  useful 
in  formulating  a  partial  differential  equation  for  the  problem.  The 
partial  differential  equation  is  obtained  by  substituting  Darcy's  law 
in  the  continuity  equation.  If  the  porous  media  is  homogeneous  and 
isotropic  and  the  flow  is  assumed  to  be  steady,  the  equation  reduces 
to  the  well-known  Laplace  equation. 

Before  attempting  to  apply  Darcy's  law  to  any  porous  media 
flow  problem,  it  is  essential  to  know  the  limitation  of  the  law. 
Darcy's  law  is  subject  to  the  following  limitations:  (1)  The  fluid  is 
homogeneous.  (2)  There  is  no  interaction  of  fluid  and  the  porous 
media.  (3)  The  flow  rates  are  relatively  small  so  that  the  inertia 
effects  are  negligible.  The  presence  of  clay  in  a  porous  media  which 
interacts  with  the  fluid  could  well  reduce  the  validity  of  the  Dareian 
proportionality.  Good  agreement  with  Darcy's  equation  occurs  for 


seepage  flow  in  nonswelling  silt-sized  particles  down  to  the  0.002 
mm  size,  generally  considered  the  upper  limit  of  the  clay  range 
(Olsen,  1965  and  1966).  However,  should  the  porous  medium 
contain  particles  within  the  clay  range,  greater-than-proportional 
response  of  flow  velocity  to  hydraulic  gradient  may  be  observed. 
Indications  suggest  the  deviation  from  linearity  is  mainly  due  to 
non-Newtonian  liquid  behavior  caused  by  clay  water  interaction 
(Swartzendruber,  1962a,  b,  1963).  Darcy's  law  is  also  valid  only  for 
low  rates  of  flow.  The  limiting  Reynolds  number  using  the  effective 
diameter  of  the  soil,  d10  as  the  length  parameter  (Re  =  Vd  10/v  in 
which  V  is  the  flow  rate  per  unit  area,  d  10is  the  effective  diameter 
of  the  opening  of  the  sieve,  through  which  10  percent  of  the  soil 
sample  by  weight  passes,  and  v  is  the  kinematic  viscosity  of  water) 
has  been  found  within  the  range  of  3  to  10.  The  limiting  Reynolds 
number  is  approximately  unity  when  using  the  average  grain 
diameter,  defined  by  d  =  31  I!  ns  d  SJ  /  Sn  s  in  which  ds  is  the 
arithmetic  mean  of  the  opening  in  any  two  consecutive  sieves  of  the 
Taylor  or  U.S.  Standard  series,  and  ns  is  the  number  of  grains  of 
diameter  ds,  as  found  by  a  sieve  analysis. 

Many  saturated  porous  media  flow  problems  have  been  solved 
using  Darcy's  law  and  the  equation  of  continuity.  Among  these  are  a 
number  of  analytic  solutions  to  a  variety  of  porous  media  flow 
problems.  (See  Muskat,  1937;  Harr,  1962;  Polubarinova-Kochina, 
1962.)  With  the  development  of  high-speed  digital  computers,  more 
complicated  problems  can  be  solved  by  numerical  methods.  A 
number  of  approximate  solutions  have  been  obtained.  (See  Shaw 
and  Southwell,  1941;  Fayers  and  Sheldon,  1962;  Freeze  and 
Witherspoon,  1966;  Zienkiewicz,  Mayer,  and  Cheung,  1966; 
Terzidis,  1968;  Jeppson,  1969a;  Taylor  and  Luthin,  1969.) 

The  problem  of  steady  infiltration  of  rainfall  into  the  surface 
of  a  watershed  is  difficult  to  solve.  The  flow  is  partially  saturated, 
and  is  effected  by  capillary  tensions  at  the  air-water  interface  as  well 
as  by  reduced  effective  soil  voids.  Modification  of  Darcy's  law  to 
apply  to  partially  saturated  flow  in  porous  media  has  been 
accomplished  by  past  research. 

An  extension  of  Darcy's  law  to  partially  saturated  flow  in  soils 
was  made  by  Buckingham  (1907)  who  studied  the  capillary  flow  of 
water  through  a  soil  with  the  analogy  of  thermal  and  electrical 
conduction.  Gardner  (1936)  and  Gardner  and  Gardner  (1950) 
suggested  that  Darcy's  equation  be  modified  for  partially  saturated 
flow  by  including  a  function  of  the  moisture  content,  f,  giving  V=  - 
Kf  grad  h.  Hall  (1956)  showed,  by  analytical  derivation,  that  the 
general  form  of  Darcy's  equation  for  anisotropic  porous  media  is 
valid  for  all  saturated  and  partially  saturated  porous  media  systems 
if  the  inertia  forces  of  the  system  are  negligible,  the  liquid  films  are 
continuous,  and  a  volume  element  of  the  fluid-porous  medium 
system  can  be  selected  which  is  small  compared  to  the  gross 
dimension  of  the  system,  yet  large  enough  that  the  surface  area  of 
the  matrix  therein  can  be  considered  to  be  uniformly  distributed 
throughout  the  volume  element.  He  also  stated  that  the  perme- 
ability K  is  a  scalar  quantity  whose  magnitude  may  vary  with  time, 
location,  and  moisture  content.  The  evaluation  of  this  quantity  by 
analytical  means  is  difficult  since  in  partially  saturated  porous  media 
flow,  capillary  forces  are  present  at  each  air-water 'interface  in  the 


interior  of  the  flow  system  and  at  the  boundary,  thus  permeability 
must  be  determined  experimentally.  The  functional  relationship 
between  saturation,  capillary  pressure,  and  the  permeability  of  a  soil 
are  important  in  obtaining  solutions  to  boundary  value  problems 
describing  the  steady  flow  of  water  through  a  partially  saturated 
porous  medium.  Gardner  (1958)  studied  available  data  and  sug- 
gested an  equation  of  the  form  K  =  a'  /(pcT  +  b1 )  where  K  is 
permeability,  pc  is  capillary  pressure  and  a',  b1,  and  r  are 
parameters  depending  on  the  liquid,  the  soil,  and  the  capillary 
pressure  history.  This  equation  has  been  used  to  solve  partially 
saturated  drainage  and  subirrigation  problems  (Sewell  and 
Schilfgaarde,  1963).  Further  investigations  made  by  Scott  and 
Corey  (1961)  and  Brooks  and  Corey  (1964,  1966)  suggested  an 
equation  of  the  form  Kr  =  (pblpc)^  for  pc  >  pb  in  which  Kr  is 
relative  permeability  defined  as  the  ratio  of  permeability  at  any 
given  saturation  to  the  permeability  at  complete  saturation,  pc  =  -  p 
is  the  capillary  pressure,  and  pb  and  r\  are  parameters  depending 
upon  the  liquid,  the  soil,  and  the  capillary  pressure  history  of  the 
system.  The  equation  ignores  the  zone  of  relatively  constant 
permeability  for  the  range  of  capillary  pressure  less  than  pb.  King 
(1965)  noted  that  Gardner's  equation  is  dimensionally  inconsistent 
and  proposed  a  dimensionless  form  of  the  equation 


(S' 


for    p  <   0 


(1) 


+    b 


in  which  K"7  is  a  dimensionless  permeability  parameter  equal  to  the 
relative  permeability  divided  by  b,  and  pb  (a  positive  quantity)  is  a 
parameter  having  the  same  dimension  as  capillary  pressure,  and  n 
and  b  are  dimensionless  parameters.  Equation  1  fits  imbibition  and 
desaturation  experimental  data  quite  well  (King,  1965).  The 
magnitudes  of  the  parameters  in  Eq.  1  are  related  to  hydraulic 
properties  of  soils.  Furthermore,  Eq.  1  can  readily  be  used  for 
solutions  of  steady-state  porous  media  flow  problems.  Also  Eq.  1  is 
algebraic  and  differentiable  and  consequently  its  use  simplifies  the 
computer  programming  for  numerical  solutions  over  that  needed  for 
a  non-algebraic  equation.  The  equation  is  also  simple,  can  be  solved 
repeatedly  without  excessive  computer  execution  time,  and  does 
not  require  core  storage  as  would  be  the  case  if  actual  data  were 
stored.  Hence  solutions  of  larger  problems  such  as  the  general 
three-dimensional  infiltration  problems  are  feasible. 

Efficient  watershed  management  requires  accurate  information 
on  the  infiltration  rate  at  which  different  soils  will  take  water  under 
different  conditions.  Infiltrometers  are  often  used  to  determine  the 
final  infiltration  rate  of  a  soil  on  small  watersheds  or  on  experi- 
mental or  sample  areas  within  large  watersheds.  However,  these 
measurements  are  subject  to  an  unknown  effect  due  to  the  lateral 
movement  of  water  in  soils.  Musgrave  (1942)  pointed  out  the 
importance  of  lateral  movement  of  water  which  occurs  to  some 
extent  with  any  method  of  applying  water  to  a  small  plot  by  rainfall 
simulators  of  various  kinds  and  sizes  and  with  different  sizes  of 
concentric  rings.  Duley  and  Domingo  (1943),  working  with  spray 
applications,  found  that  the  mean  total  intake  of  water  on  isolated 
small  plots  (16  by  72  inches)  having  no  prewetted  border  protection 
was  75  percent  higher  than  for  the  large  plots  protected  with  a 
buffer  strip  and  a  saturated  belt  of  soil,  and  92  percent  greater  than 
for  the  small  enclosed  plots.  The  results  indicated  that  lateral 
movement  of  moisture  had  allowed  these  small  plots  to  take  in  more 
water  than  would  be  possible  if  it  were  raining  over  the  entire 
surface  of  a  large  area.  Marshall  and  Stirk  (1950)  used  both 
buffered  and  unbuffered  rings  ranging  from  1  to  10  feet  in  diameter 
to  study  the  effect  of  size  of  plots  on  the  lateral  movement  of  water 
and  showed  that  the  steady  minimum  infiltration  rate  of  a  given  soil 
decreased    with    increasing    size   of  plot.    Wooding   (1968),   upon 


assuming  that  permeability  can  be  defined  as  an  exponential  of  the 
pressure,  linearized  the  non-linear  differential  equation,  and  ob- 
tained solutions  for  steady  infiltration  from  a  shallow,  circular, 
inundated  area  on  the  horizontal  surface  of  a  semi-infinite  porous 
medium,  and  showed  varying  degrees  of  lateral  movement  of 
moisture  for  varying  soils.  Jeppson  and  Nelson  (1970),  using  Eq.  1 
to  describe  hydraulic  properties  of  soils  for  partially  saturated 
seepage  from  canals,  obtained  solutions  and  found  that  a  more 
pronounced  spreading  effect  of  moisture  movement  can  be  detected 
for  the  partially  saturated  solutions  than  for  solutions  based  on  the 
saturated  flow  assumption. 

Investigations  on  the  spreading  effect  and  infiltration  charac- 
teristics through  the  application  of  partially  saturated  flow  theory 
to  infiltration  is  important  to  the  development  of  watershed 
hydrology.  To  investigate  spreading  characteristics  is  one  of  the 
objectives  of  this  study. 

The  problem  of  infiltration  of  rainfall  on  the  surface  of  a 
watershed  can  be  formulated  mathematically  by  non-linear  partial 
differential  equations.  The  general  three-dimensional  unsteady 
problem  thus  formulated  is  difficult  to  solve.  The  axisymmetric  case 
has  been  solved  by  Jeppson  (1970c).  The  problem  of  steady-state 
axisymmetric  infiltration  resulting  from  moisture  applied  at  a 
horizontal  surface  has  been  selected  for  this  study  to  determine  the 
influence  of  various  soil  properties  (i.e.  the  parameters  in  Eq.  1 
which  describe  a  wide  variety  of  soils  and  their  hydraulic  properties) 
on  the  flow  patterns  as  well  as  the  magnitude  of  the  "spreading 
effect"  due  to  the  lateral  movement  of  soil  moisture. 

The  steady-state  infiltration  rate  which  represents  the  final 
long-term  intake  rate  at  which  water  will  be  absorbed  by  the  plot 
area  is  important  in  watershed  management.  The  information 
obtained  from  the  steady-state  solution,  will  give  guidance  in 
estimating  effects  of  long  duration  rainfalls  on  the  watershed.  In 
addition  to  giving  information  concerning  the  relative  importance  of 
the  "spreading  effect,"  the  steady-state  solutions  will  provide  a 
check  on  the  validity  of  unsteady-state  solutions. 

The  partial  differential  equations,  which  result  from  substitut- 
ing Eq.  1  into  the  steady  state  continuity  equation  as  well  as  the 
equivalent  equation  using  the  stream  function  as  the  dependent 
variable,  are  of  elliptic  type.  With  proper  boundary  conditions, 
axisymmetric  flow  from  an  infiltrometer  becomes  a  mathematical 
boundary  value  problem.  Reisenauer  (1963),  using  the  methods 
outlined  by  Nelson  ( 1962)  has  obtained  a  numerical  solution  to  the 
problem  of  steady-state  partially  saturated  axisymmetric  flow  from 
a  circular  pond  by  a  formulation  in  the  physical  plane.  An  inverse 
formulation  (see  Jeppson,  1966,  196^7,  1968a,  1968b,  1969a, 
1969b;  Jeppson  and  Nelson,  1970)  has  been  used  in  this  study  to 
construct  the  mathematical  problem.  This  inverse  formulation 
considers  the  coordinates  in  the  radial  and  axial  directions  (r  and  z 
respectively)  as  the  dependent  variables  and  the  potential  function 
O  (proportional  to  the  hydraulic  head)  and  the  stream  function  ^ 
as  the  independent  variables  (i.e.  the  problem  is  solved  for  r  and  z  in 
the  $j;  plane).  Finite  difference  techniques  are  used  in  solving  the 
problem  by  an  iterative  procedure  and  a  computer  program  has  been 
written.  The  solution  results  consist  of  the  r  and  z  coordinates  at  the 
intersection  of  each  streamline  with  each  potential  line.  These  values 
can  readily  be  used  to  construct  the  flownet  with  a  computer  driven 
plotter.  The  resulting  flownets  show  characteristics  of  the  flow 
patterns  at  a  glance.  In  addition,  the  capillary  tension  can  be 
obtained  at  each  point  from  the  values  of  z  and  the  potential  $ 
from  the  solution,  and  the  relative  permeability  at  each  point  can  be 
obtained  by  substituting  capillary  tension  at  each  point,  and  the  soil 
parameters  specified  to  obtain  the  solution,  into  Eq.  1. 


Saturation    is  another  item   of  interest.  The   functional  rela-  cal  approximation  methods  have  been  used  to  obtain  values  for  Sc 

tionship    between    relative    permeability    and    saturation    can    be  corresponding  to  Kr. 
obtained    analytically,    if  a   suitable   algebraic   expression   for  the 

effective  saturation  Se  =  f(pc)  can  be  defined,  which  involves  the  In  the  following  sections,  the  mathematical  formulation  and 

soil    parameters   in    Eq.    1.   While   this  algebraic  expression  is  not  the  procedure  followed  in  obtaining  the  inverse  finite  difference 

available,    the    first  order  ordinary  differential  equation  obtained  operators  for  the  problem  and  the  numerical  techniques  involved  in 

from  Burdine's  equation  (Burdine,   1953)  for  relative  permeability  obtaining  solutions  for  different  soil  types  are  discussed  in  detail, 

can  be  solved  for  the  relationship  Kr    vs.  Se.  The  solution  of  the  Solution  results  are  presented  in  flownet  form,  and  effects  due  to 

differential  equation  cannot  be  obtained  analytically,  thus,  numeri-  different  soil  parameters  are  given. 


DIFFERENTIAL  EQUATIONS 


The  general  form  of  the  Darcy's  equation  is  assumed  to  apply 
for  partially  saturated  flow 


K  grad 


(2) 


The  differential  equation  for  either  the  potential  function  or  the 
stream  function  can  be  obtained  from  Eqs.  5  and  6.  Upon 
differentiating  Eq.  5  with  respect  to  r  and  Eq.  6  with  respect  to  z 
and  combining  the  resulting  equations  to  eliminate  the  common 
term  S^/S  rSz  =  S2i|i  /S  zSr  gives 


in  which  v  is  the  velocity  vector,  K(p(r,z))  is  the  permeability 
coefficient  of  the  soil  with  dimension  of  velocity,  and  0  is  the 
potential  energy  per  unit  pound  of  water  with  dimension  of  length 
and  is  defined  by 
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The  general  form  of  Eq.  8,  not  associated  to  a  particular  coordinate 
system,  can  be  obtained  by  substituting  Darcy's  Eq.  2  into  the 
steady-state,  continuity  equation. 


in  which  p  is  the  pressure  of  water,  p  is  the  mass  density  of  water,  g 
is  the  acceleration  of  gravity,  and  z  is  the  vertical  coordinate.  The 
effects  of  the  soil  and  water  are  incorporated  in  the  coefficient  of 
permeability 


(4) 


in  which  k(p(r,z))  is  the  effective  intrinsic  permeability  with 
dimensions  of  length  squared,  and  jjl  is  the  absolute  viscosity  of 
water.  The  velocity  components  also  result  from  the  definition  of 
Stokes'  stream  function.  Thus  the  components  of  velocity  in  the 
radial  and  axial  directions  are  given  by 
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V  •  (K  grad  0)     =     K  V    0   +  grad  K.grad  0     =     0 


(9) 


The  equation  for  the  stream  function  which  results  from 
differentiating  Eq.  5  with  respect  to  z  and  Eq.  6  with  respect  to  r 
and  combining  the  resulting  equations  is 
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Equations  8  and  10  apply  in  both  the  saturated  and  the 
partially  saturated  regions  of  flow.  In  saturated  regions,  since  K  is 
constant  for  homogeneous  soil,  the  terms  including  derivatives  of  K 
vanish.  In  partially  saturated  regions,  the  derivatives  with  respect  to 
K  must  be  evaluated  from  the  relationship  between  permeability 
and  capillary  pressure.  The  equation  used  for  this  relationship  is 
described  later. 


and 


Instead  of  solving  Eqs.  8  and  10  on  the  physical  plane  as 
mentioned  earlier,  the  problem  in  this  study  is  solved  on  the  Oij/ 
plane  (a  rectangular  region). 
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In  Eqs.  5  and  6,  it  is  assumed  that  the  flow  is  axisymmetric  and  the 
soil  is  isotropic,  i.e.  the  permeability  is  the  same  in  the  radial  and 
axial  directions.  For  anisotropic  soils,  with  the  permeability  in  the 
radial  direction  equal  to  Kr  and  that  in  the  vertical  direction  Kz, 
Eqs.  5  and  6  apply  if  z  is  replaced  by  zt  and  K  by  Kt  from  the 
following  two  equations  (Muskat,  1937). 
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Inverse  formulation 

The  inverse  partial  differential  equations  which  describe  the 
flow  in  terms  of  r(4> ,  jj)  and  z(cj>,ij;)  can  be  derived  from  the 
implicit-function  theory  (Taylor,  1955)  following  the  approach  used 
by  Jeppson  ( 1968b).  A  pair  of  basic  equations  similar  to  Eqs.  5  and 
6  can  be  obtained; 
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These  equations  are  important  for  developing  the  governing  partial 
differential  equations  in  the  Ody  plane.  Upon  differentiating  Eq.  1 1 
with  respect  to  0  and  Eq.  1  2  with  respect  to  ty,  and  combining  the 
results  to  eliminate  the  common  term  d2  z/d^d^i  =  B2  z/dij;d<t> 
gives 


in  which  Kr   is  the  permeability  parameter  defined  earlier  and  given 
by 
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Differentiating  Eq.  16  with  respect  to$  and  4*  gives  respectively 
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Also  differentiating  Eq.    1 1    with  respect  to  i|j    and  Eq.    1 2  with 
respect  to  $  and  combining  the  results  gives 
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Since  both  dependent  variables  r(0,i|i)  and  z($,ijj)  appear  in  Eqs. 
13  and  14,  these  two  equations  must  be  solved  simultaneously.  In 
partially  saturated  regions  K  is  a  function  of  pressure,  and 
consequently,  also  a  function  of  z.  The  methods  used  for  evaluating 
K  and  its  derivatives  will  be  discussed  in  the  following  section. 

Evaluation  of  the  derivatives 
of  the  permeability 

Equation  1  proposed  by  King  (1965)  which  describes  the 
relationship  between  capillary  pressure  and  permeability  at  any 
given  saturation  for  imbibition  as  well  as  drainage  was  adopted  for 
this  study.  For  convenience  Eq.  1  is  rewritten 


(14)  _ 


K 


K    = 


(1) 


+  b 


in  which  KQ/b  is  the  saturated  permeability,  p(r,z)  is  the  water 
pressure  and  is  negative  for  partially  saturated  flow,  pb  is  the 
bubbling  pressure  and  r\  is  the  pore  size  distribution  index  defined 
as  the  negative  slope  of  the  straight-line  portion  of  the  log-log  plot 
of  Kr  versus  (-p/pj,)  and  b  is  a  dimensionless  constant.  The 
magnitude  of  the  parameters  in  Eq.  1  depends  on  the  characteristics 
of  the  soil  and  the  capillary  pressure  history  of  the  soil  water 
system. 


_        1 

dK 

r              v 

PfJ  -^'  1 

dijj 

}$  e^  ♦  b' 

2 

,,-2  IPs  ^      ,      ^Tl-1    dz 


\Pb/ 
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(19a) 


(19b) 


Eqs.  15  through  19b  apply  only  in  those  regions  of  flow  which  are 
partially  saturated  (i.e.  when  p  <  0).  Whenever  the  soils  become 
saturated,  the  coefficient  of  permeability  remains  constant  and 
equals  K0/b. 

Scaling  of  variables 

Since  permeability  is  generally  a  relatively  small  quantity,  Eqs. 
13  and  14  can  be  made  more  amenable  to  solution  by  finite 
difference  methods  by  making  the  coefficients  of  each  term  closer 
to  unity.  To  accomplish  this,  the  total  potential  energy  per  pound 
of  water  O  has  been  scaled  by 


>      =    K    <2> 
t  o 


(20) 


Using  the  definition  of  the  total  energy  per  pound  of  water 
from  Eq.  3 


From  Eq.  20,  the  following  relationships  can  be  obtained 


Pg(z  -*) (15) 


Substituting  Eq.  1  5  into  Eq.  1  gives 
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.(22) 
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•  (23)  Equations  26,  27  and/or  28,  29  are  final  equations  for  governing  the 
partially  saturated-saturated  flow  system  in  the  transformed  $,  ty 
plane. 


bt 


=     K 


.  (24)         Boundary  conditions 

Equations  26  and  27  are  elliptic  types.  (See,  e.g.,  Courant  and 
Hilbert,    1962,    for    definition    of   elliptic    type    equation.)    Con- 
sequently, the  problem  being  studied  is  a  boundary  value  problem. 
(25)  For  a  well  posed  boundary  value  problem,   boundary  conditions 

must  be  established  on  all  boundaries  enclosing  the  region  of  flow. 


Upon  substituting  from  Eqs.  21  through  25.  Eq.  13  becomes 


—        /— 2      2     S2r        a2r 
K      r    K       r        — 7   +  5 

r      \  r  do,2      exs2 
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The  problem  which  is  described  herein  is  that  of  steady  state 
axisymmetric  infiltration  from  a  circular  rainfall  simulator  through 
homogeneous  soil  to  a  water  table  at  a  known  depth  (see  Fig.  1). 
Since  the  rainfall  is  applied  over  a  circular  area  of  the  soil  surface, 
the  region  of  infiltration  below  the  circular  area  is  symmetric  about 
the  vertical  centerline.  The  boundary  value  problem  can  be  set  up 
for  only  one-half  of  any  vertical  plane  containing  the  axis  of 
symmetry  such  as  the  vertical  rectangular  area  shown  in  Fig.  1. 


„    „   3k     -,       Sk     -. 

—2     2  r   o_r         r    or 

|Kr    f        04;     04,  "  SO       0$ 


+  K 


Kr  r  \bl)  "  (a*  ) 

L  * 


.(26) 


The  boundary  conditions  are  shown  in  Fig.  2  in  which  the 
equivalent  flow  region  to  that  shown  in  Fig.  1  has  been  plotted  on 
the  Ot  4*  plane.  These  boundary  conditions  are  used  to  establish 
values  of  variable  r  along  the  water  table  and  the  soil  surface 
including  the  circular  area  of  surface  infiltration  and  the  values  of 
the  variable  z  along  the  axis  of  symmetry.  These  conditions  are 
obtained  as  follows  (the  number  notations  for  the  boundaries  are 
shown  on  Figs.  1  and  2). 


Likewise,  Eq.  14  becomes 


, — 2     2    o    z    ,    o    z 
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r       \    r  Sty2       oi2 


—  2     2       r  dz         r    oz 
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Axis  of  symmetry  Q)  -  Q).  The  value  of  r  is  zero  along  the  axis 
of  symmetry  and  since  equipotential  lines  are  horizontal  in  the 
physical  plane  at  the  axis  of  symmetry,  the  boundary  conditions 
along  the  axis  of  symmetry  are  r  =  0  and  o  z/o^    =  0. 

Special  consideration  must  be  given  in  obtaining  the  solution 
near  the  axis  of  symmetry,  since  as  the  radius  r  approaches  zero,  the 
inverse  transformation  used  to  obtain  Eqs.  1 1  and  1  2  is  not  valid. 
Failure  of  the  transformation  at  r  =  0  can  be  checked  by  noting  that 
the  value  of  the  Jacobian  J  for  the  transformation,  which  is  given  by 


.(27) 


—  2    o_z     6z 
+  2Kr      04;    0$ 


In  the  case  of  saturated  soils,  the  permeability  parameter  Kr 
remains  constant  and  equals  1/b.  The  derivatives  of  Kr  with  respect 
to  both  $,  and  6  are  consequently  zero  whenever  p  5  0  and  Eqs. 
26  and  27  reduce  respectively  to  the  inverse  equation  for  axi- 
symmetric seepage  through  saturated  homogeneous  porous  media 
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approaches  zero  as  r  approaches  zero,  with  dz/c^   =  0,  i.e. 
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and 


2  a  r      a  z     a_z  az 
r   34,2+  »,J+  a**, 


=    0 


Equation  31  shows  that  a  line  singularity  exists  along  the  axis  of 
symmetry  and  the  derivative   Sr/Sijj    is  discontinuous  at  the  axis. 
.(29)  -T-^g  manner  jn  which  this  line  singularity  has  been  handled  in  the 

finite  difference  solution  is  described  in  a  later  section. 
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Fig.  1 .      Axisymmetric  infiltration  through  partially  saturated  porous  media  from  surface  applica- 
tion of  moisture  over  a  circular  area  to  a  water  table. 
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Fig.  2.      Formulation  of  the  boundary   value  problem  for  steady  state  axisymmetric  infiltration 
from  a  source  circle  through  a  homogeneous  porous  media  to  a  water  table. 


Source  surface  Q)  -  Q).  The  source  surface  is  assumed  to  be 
horizontal  and  consequently  at  a  constant  height  D  above  the  water 
table,  the  condition  for  z  along  this  boundary  is 


z    =    D  , 


.(32) 


provided  that  the  origin  of  the  coordinates  is  selected  at  the  axis  of 
symmetry  on  the  water  table.  A  constant  pressure  head  is  also 
specified  over  the  source  surface  resulting  in  a  constant  value  for  K. 
Consequently  the  potential  energy  per  pound  of  fluid  at  the  surface 
is  also  constant,  i.e. 


NH 


H 


.(33) 


The  value  of  H  is  read  as  an  input  parameter  in  the  computer 
program.  If  the  value  of  H  is  selected  to  be  greater  than  D,  ponded 
water  with  depth  H-D  exists  over  the  source  area,  and  the  problem 
becomes  a  combined  saturated-unsaturated  flow  case.  Also  since  the 
streamlines  leave  the  source  surface  vertically,  the  boundary 
condition  for  r(  0,  ,\\i )  can  be  expressed  by 


c)<J> 
t 


.(34) 


perpendicularly,  the  condition  for  r  is  Sr/c)4/  =  0.  The  boundary 
condition  for  r  can  also  be  obtained,  utilizing  the  solution  forz,  by 
integrating  Eq.  II  along  the  boundary  (i.e.  the  last  ijj  =  constant 
line),  giving 


--.I 


k"    r    U    d* 
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.(36) 


The  subscript  by  the  integral  indicates  the  integration  is  to  be 
carried  out  along  a  streamline.  This  implicit  equation  for  r  is  solved 
by  finite  differences. 

Vertical  outer  boundary  Q)  -  Q).  This  vertical  boundary  is 
taken  at  a  great  enough  distance  from  the  axis  of  symmetry  so  that 
no  moisture  movement  will  occur  along  the  boundary.  The 
boundary  is  considered  as  a  streamline  and  equipotential  lines  are 
essentially  horizontal  at  this  boundary.  Under  this  assumption  the 
condition  becomes  bz/b^i  =  0.  An  alternative  approach  may  be 
used  to  obtain  values  of  z  on  and  adjacent  to  part  of  this  boundary. 
This  approach  has  not  been  used  in  this  study,  but  reduces  the 
partial  differential  equation  to  an  ordinary  differential  equation  by 
noting  that  Kr  is  very  small  along  most  of  this  boundary. 
Consequently  terms  in  Eq.  27  which  contain  K^  2  or  Kr  3  can  be 
omitted.  Equation  27  then  becomes 


Values  of  r(  J,.^)  along  the  source  surface  boundary  may  be 
established,  alternatively,  by  integrating  Eq.  12  (Jeppson,  1968b). 
After  transforming  $  to   <£,  and  integrating.  Eq.  1  2  yields 
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The  subscript  of  the  integral  indicates  the  integration  is  to  be  carried 
out  along  an  equipotential  line.  i.e.  along  the  boundary.  This 
implicit  equation  for  r  (notice  that  r  appears  on  both  sides  of  the 
equal  sign)  can  be  solved  by  finite  differences. 

Ground  surface  ©  -  @.  The  ground  surface  is  horizontal  and 
at  a  constant  height  D  above  the  water  table,  hence,  the  condition 
for  z  along  this  boundary  is  z  =  D.  Since  the  ground  surface  is  a 
streamline,    with   each   equipotential   line   entering   the   boundary 


Values   of  z   along   the   outer  vertical   boundary  come  from   the 
solution  of  Eq.  37. 

Water  table  Q)  -  Q).  The  water  table  is  specified  horizontal. 
The  boundary  condition  for  z  along  this  boundary  is  z  =  0.  Since  the 
streamlines  enter  the  water  table  vertically,  the  boundary  condition 
for  r(  $t  \\i  )  can  be  expressed  by  Sr/c)  4*t  =  0-  The  alternative  way 
to  establish  the  values  of  r($  t  ^ )  along  the  boundary  by  integration 
as  mentioned  in  the  last  section  has  not  been  used,  because  the  rapid 
variations  in  Kr  in  the  vicinity  of  the  boundary  have  caused 
difficulty  in  obtaining  a  solution  of  r  on  the  boundary.  Instead,  a 
parabolic  extrapolation  has  been  used  to  meet  the  normal  condition 
Sr/S$t    =  0  and  will  be  described  later. 


INFLUENCE  OF  PARAMETERS  IN  PERMEABILITY-CAPILLARY 
PRESSURE  RELATIONSHIP 


In  partially  saturated  flow  in  porous  media,  capillary  forces  are 
present  at  each  air-water  interface  in  the  interior  of  the  flow  system. 

Several  empirical  equations  have  been  proposed  for  representing  the 
permeability  as  a  function  of  capillary  pressure  in  terms  of  some  soil 
parameters.  The  manner  in  which  the  permeability  varies  with  the 
parameters  in  two  of  these  equations  is  described  in  this  section. 


Scott  and  Corey  (1961)  proposed  the  equations 
for    p     >  p        


K 


and 
K 


PbP 


for     p^  <  p. 


.(38) 


in  which  Kr  is  the  relative  permeability  (the  ratio  of  permeability  at 
any  given  saturation  to  the  permeability  at  complete  saturation).  pc 
is  the  capillary  pressure  (the  absolute  value  of  the  pressure 
difference  between  air  and  water  at  the  air-water  interface  in  the 
pores  of  soil),  pb  is  bubbling  pressure  and  was  found  by  Brooks  and 
Corey  (1964)  to  be  closely  related  to  the  largest  pores  forming  a 
continuous  network  within  a  porous  medium,  and  r\  is  the  pore-size 
distribution  index  which  is  expressed  by  the  negative  slope  of  the 
straight  line  portion  of  the  log-log  plot  of  Kr  versus  pc .  Equation  38 
defines  a  straight  line  on  log-log  paper  having  an  intercept  at  Kr  = 
1.0.  This  equation  assumes  the  permeability  is  constant  for  the 
range  of  capillary  pressure  less  than  pb  (Brooks  and  Corey,  1964, 
and  King,  1965). 

King  has  defined  this  zone  of  relatively  constant  permeability 
as  a  "plateau"  which  may  exist  for  both  drainage  and  imbibition 
(King,  1965).  King  has  also  proposed  Eq.  1  as  a  dimensionless  form 
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Fig.  3.      Permeability  parameter,  Kt,  as  a  function  of  capillary  pressure  head  for  different  pb/7 
values,  with  r|  =  2.50  and  b  =  1.0. 


of  Gardner's  equation  as  mentioned  earlier.  Equation  1  and  not  Eq. 
38  was  used  to  represent  the  permeability-capillary  pressure 
relationship  in  this  study  and  is  discussed  in  detail  in  this  section. 


The  functional  relationship  between  Kr  and  capillary  pressure 
head  obtained  from  Eq.  1  is  shown  in  Fig.  3  for  different  values  of 
pb/7  with  tj  =  2.50  and  b  =  1.0.  The  five  curves  on  the  graph  are 
nearly  parallel  for  larger  values  of  pcly.  The  lengths  of  plateaus  are 
different  depending  on  pb/Y  and  to  a  lesser  extent  on  b.  The  length 
of  the  plateau  decreases  as  the  value  of  bubbling  pressure  decreases. 
Curves  for  TT7  as  a  function  of  capillar,  pressure  are  shown  in  Fig.  4 
for  different  values  of  b  with  pb/y  =  3.50'  and  n  =  2.50.  The 
slopes  of  the  straight  line  portion  of  the  curves  are  the  same. 
Furthermore,  the  capillary  pressures  obtained  by  extrapolation  of 
the  straight  line  portion  of  the  curves  at  the  top  of  the  graph  where 
Kr  =  1.0  are  the  same.  The  main  difference  in  these  curves  is  the 
height  of  plateau,  which  decreases  with  the  increasing  value  of  b. 
Five  K7   versus  pc    curves  obtained  from  Eq.  1  are  shown  in  Fig.  5 


for  different  n-values  with  pb  /7  =  3.50'  and  b  =  1.0.  The  plateaus 
of  these  curves  have  the  same  heights,  and  the  capillary  pressures 
obtained  by  extrapolation  of  the  straight  line  portion  of  the  curves 
to  the  value  of  Kr  =1.0  are  also  the  same,  and  equal  to  the 
bubbling  pressure.  On  the  graph,  all  the  curves  have  passed  through 
a  common  point  at  pc  =  pb  and  Kr  =  l/(l+b).  For  pc  >  pb  ,  the 
permeability  parameter  (Kr  )  increases  as  the  value  of  "H  decreases 
for  constant  capillary  pressure  heads  greater  than  pb  jy  and  Kr 
decreases  as  the  value  of  r\  decreases  for  pc  <  pb  .  Also  the  radius 
of  curvature  of  the  curved  portion  of  the  curve  increases  as  the  value 
of  1  decreases. 

The  function  Kr  vs.  p/7  will  approach  a  step  function  as  the 
value  of  i  becomes  very  large. 

The  properties  of  Eq.  1  depicted  in  Figs.  3,  4,  and  5  have 
influences  on  the  flow  patterns  resulting  from  specification  of 
different  values  of  the  parameters  b,  pb  and  r|.  These  influences  are 
discussed  in  later  sections. 
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Fig.  4.     Permeability  parameter,  Kr ,  as  a  function  of  capillary  pressure  head  for  different  b  values, 
with  pb/7  =  3.50'  and  n  =  2.50. 
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Fig.  5.     Permeability  parameter,  Kr ,  as  a  function  of  capillary  pressure  head  for  different  *1  values, 
with  pb/-y  =  3.50'  and  b  =  1 .0. 
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RELATIVE  PERMEABILITY  AS  A  FUNCTION  OF  SATURATION 


If  air  and  water  are  the  nonwetting  and  wetting  fluids. 
respectively,  experimental  data  of  effective  saturation  vs.  capillars 
pressure  plots  close  to  a  straight  line  on  log-log  graph  paper  for 
capillary  pressure.  pc  .  greater  than  bubbling  pressure.  pb  .  A  simple 
equation  defining  this  relationship  for  the  drainage  cycle  was  given 
by  Brooks  and  Corey  ( 1 964) 


The  necessary  Burdine  integration  has  been  carried  out  for  Eq. 
39  (or  Eq.  41  with  A  =  0)  by  Brooks  and  Corey  (1964)  to  obtain 
the  approximation  equation  for  relative  permeability,  i.e.  Eq.  38. 
And  the  relationslup  for  Kr  vs.  Se  becomes  a  function  of  the 
pore-size  distribution  index,  n,,  only,  and  the  following  approxima- 
tion equation  for  Kr   vs.  Se   is  obtained 


for    P_>Pb 


.(39) 


3X+2 
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where  A    =  (r|-2)/3   is  a  number  closely   related  to  the  pore-size 
distribution  index,  r|.  and  the  effective  saturation  is  defined  by 


1   -  S 


.(40) 


When  A  is  not  zero,  the  limit  for  Se  approaches  the  constant 
value  1/A  as  pc  approaches  zero.  Generally,  letting  A  take  on  a 
value  gives  a  better  fit  to  experimental  data  than  if  A  =  0.  With  a 
non  zero  A  in  Eq.  41  there  is  no  direct  method  such  as  integrating 
the  Burdine  equation  for  evaluating  the  magnitudes  of  A  and  A  (Eq. 
42)  from  the  magnitudes  of  the  parameters  in  Eq.  1  as  is  the  case  if 
A  =  0.  Consequently,  the  Burdine  equation  has  been  solved 
numerically. 


in  which  S  is  the  saturation  and  Sr  is  the  residual  saturation.  (See 
Brooks  and  Corey,  1964,  for  the  determination  of  Sr  .)  A  more 
general  approximation  equation  proposed  by  Brutsaert  ( 1968)  is 


.(41) 
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where  A  and  A1  are  the  parameters  whose  magnitude  depends  on 
the  soil  types.  Where  A  is  zero,  the  Brooks-Corey  approximation 
results. 


If  Eq.  42  is  differentiated  with  respect  to  Se  ,  a  first  order 
differential  equation  can  be  obtained 
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According  to  the  analysis  of  Burdine  (1953),  the  relative 
permeability  can  be  expressed  as  a  function  of  Se  (Brooks  and 
Corey,  1964) 


is  a  constant  depending  on  the  relationship  of  saturation  to  capillary 
pressure.  While  Eq.  44  cannot  be  solved  analytically  the  solutions  of 
the  equation  can  be  generated  numerically  if  a  correct  starting  value 
can  be  found.  The  expression  for  Kr    is  given  by 
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An  approximate  value  of  C  can  be  obtained  by  integrating  Eq.  39 
giving 


where  Kr    is  the  relative  permeability  which  has  a  maximum  value 
of  unity  as  pc    approaches  zero. 
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(46) 


12 


With  this  approximate  value,  Eq.  44  can  be  solved  approximately 
and  thereafter  the  value  of  C  can  be  improved  by  using  the  resulting 
relationship  for  Se  and  pc  .  This  improved  value  of  C  can,  in  turn, 
again  be  used  to  obtain  a  better  solution  to  Eq.  44.  By  repeating  this 
procedure  until  the  value  of  C  does  not  change  between  consecutive 
cycles,  the  solution  is  obtained.  And  the  solution  to  Kr  vs.  Se  is 
obtained  from  the  solution  to  Eq.  44  when  Kr    =  Kr  /b. 

To  solve  Eq.  44  it  is  necessary  to  start  with  a  value  of  K  r 
which  corresponds  to  the  starting  value  of  Se  .  When  the  capillary 
pressure  approaches  infinity,  the  permeability  parameter,  Kr  ,  as 
well  as  the  effective  saturation  approaches  zero.  If,  however,  Kr  =0 
and  Sc  =  0  are  selected  as  the  starting  values,  the  first  term  on  the 
right  side  of  the  equal  sign  in  Eq.  44  becomes  indeterminate,  and 
the  generation  of  the  solution  cannot  be  carried  out.  To  prevent  this 
indeterminancy,  both  Kr  and  Se  have  been  assigned  small  values 
with  Kr  less  than  Se  .  The  magnitude  of  these  small  values  must  be 
improved  by  a  trial  and  error  procedure  so  that  Kr  =1.0  when  Se 
=  1.0.  The  solution  might  also  be  obtained  starting  with  Kr  =1.0 
and  Se  =  1.0  and  using  a  negative  increment  for  Se  .  If  this 
approach  is  used,  however,  it  is  necessary  to  modify  Eq.  44  to 
prevent  division  by  zero  because  pc  =  0  when  Se  =  1.  To  prevent 
the  integral  from  becoming  infinite  as  pc  approaches  zero,  the 
numerical  integration  must  stop  (or  start)  with  Se  less  than  unity. 
Alternatively,  a  constant,  p0  ,  might  be  added  to  pc  (Jeppson, 
1970b).  If  p0  is  added  to  pc  only  for  a  portion  of  the  solution 
where  1  >  Sc  >  S0  in  which  S0  is  a  constant  ranging  from  0.9  to 
1.0,  Eq.  44  becomes 


■r 


dS 


(P.  +  P    >' 


is  a  constant  depending  on  the  relationship  of  saturation  to  capillary 
pressure. 

A  semi-log  plot  is  shown  in  Fig.  6  for  the  permeability 
parameter,  Kr  ,  vs.  effective  saturation,  Se  ,  obtained  by  solving  Eq. 
44  numerically  for  different  q  values  with  pb  /y  =  3. 50'  and  b  = 
1.0  using  the  first  approach  mentioned  earlier,  that  is,  starting  with 
small  values  of  KT  and  Se  .  The  permeability  parameter  increases  as 
the  value  of  the  pore-size  distribution  index,  r| ,  increases  at  a  given 
effective  saturation.  In  other  words,  at  a  given  saturation,  when  the 
pore-size  distribution  index  increases,  the  resistance  of  the  soils  to 
the  flow  is  reduced  and  the  permeability  is  also  increased. 

Curves  for  the  permeability  parameter,  K  r  ,  vs.  effective 
saturation,  Se  ,  for  different  bubbling  pressure  heads,  pb  /  y,  with  q 
=  2.50  and  b  =  1.0  using  the  first  approach  were  also  obtained. 
Within  the  accuracy  of  the  numerical  methods,  no  differences  in  the 
Kr  ,  Se  relationships  for  different  values  of  pb  ly  could  be 
detected  and  consequently  it  is  concluded  that  the  effect  of 
bubbling  pressure  upon  the  Kr  vs.  Se  curves  is  very  small.  In  other 
words,  the  relationship  between  relative  permeability  or  the 
permeability  parameter  and  effective  saturation  is  quite  independent 
of  the  bubbling  pressure.  The  same  conclusion  can  be  drawn  from 
Eq.  43. 


dK  2K                     /       S 
r            r  _1_     I          e 

dS  S  C      \p     +  p 

e  e                     \    c          o 


for     1    >    S       >     S 


.(47) 


Figure  7  shows  curves  for  relative  permeability  vs.  effective 
saturation  for  different  b  values  with  pb  /y  =  3. 50' ,  and  r|  =  2.50 
(using  the  first  approach  to  obtain  the  solution  curves).  The  curves 
on  this  figure  show  that  for  Se  constant,  Kr  decreases  as  the  value 
of  b  decreases.  When  the  value  of  Kr  is  scaled  to  Kr  as  defined  by 
Eq.  45,  curves  similar  to  those  on  Fig.  7  would  almost  coincide  with 
each  other.  The  only  parameter  which  significantly  affects  the 
curves  for  Kr    vs.  Se  is  the  pore-size  distribution  index  q. 


in  which 


The  distribution  of  effective  saturation  for  each  solution  can 
be  obtained  from  the  r  and  z  coordinates  given  as  computer  output 
by  using  Eq.  16  and  the  curves  obtained  for  K  r  vs.  Se  . 
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Fig.  6.  Theoretical  curves  for  permeability  parameter,  K  ,  as  a 
function  of  effective  saturation  for  different  r\  values,  with 
pb/7  =3.50'  andb=  1.0. 


Fig.  7.  Theoretical  curves  for  permeability  parameter,  Kt,  as  a 
function  of  effective  saturation  for  different  b  values,  with 
p  /y  =  3.50'  and  r\  =  2.50. 
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FINITE  DIFFERENCES 


The  finite  difference  operator  for  Eq.  26  can  be  obtained  by 
replacing  the  derivatives  by  second  order  central  differences  and 
writing  the  result  as  a  polynomial  in  terms  of  r  at  the  grid  point  in 
question.  The  operator  for  Eq.  26  is 


fl<Vo>    =     2(Kr>o2    ro4-(Vc2<Vr4)ro3 


JR<Kr>o2(^)     ^V^^VVW 


+  7   (K    )       (r    -r    )2-2R2 
4         r  o        2      4 


;k- 


R   (r1+r3) 


Similar  notations  apply  for  z.  The  subscript  i  denotes  the  number  of 
constant  §tgrid  lines  (i.e.  0,  =  (i-l)A$t )  and  j  denotes  the  number 
of  the  constant  ij;  grid  lines  (i.e.  i\i  =  (j-1)A4j).  The  ratio  of  the 
spacing  of  streamlines  to  the  transformed  equipotential  lines  is 
denoted  by  R  =  AiJj/A  0, . 

For  axisymmetric  infiltration  through  saturated  homogeneous 
porous  media,  the  operators  (Eqs.  48  and  49)  reduce  to 

r  4  -\   (r,  +  r    )r  3 
o         2        2        4     o. 


+  LR     -8<*2"4> 


2       1 


R     (r,  +r    )r 
1         Jo 


2  r  o   K       \Pk/  o      o  1       1   J      o 


+  iR2(ri  "r3>2     =     ° 


.(50) 


+  i  R2<rrr3)2  =  ° 


.(48) 


and 


The  operator  for  Eq.  27  is 


°         2(R2+r2) 


f,  (r    .  z    ) 
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+  ^r-   (z.  +z_  -2z    )+-R(z     -z    )(z     -z    )    =    0 
13  o2  13Z4 

(49) 


(KA 


These  latter  two  operators  can  also  be  obtained  from  Eqs.  28  and 
29.  Only  Eq.  51  of  the  four  operators,  Eqs.  48  through  51,  can  be 
solved  explicitly  for  the  variable  z  at  the  grid  point  (r0  ,  zQ).  To 
evaluate  r  and  z  at  the  grid  point  in  question  in  the  other  operators, 
it  is  necessary  to  solve  the  operators  by  implicit  methods.  Jeppson 
(1968b,  1968c)  used  the  Newton-Raphson  iterative  method  for 
solving  implicit  functions  for  inner  iteration  and  successive  over 
relaxation  for  the  outer  iteration  to  settle  the  values  of  the  variables 
throughout  the  flow  field.  The  same  approach  has  been  used  to 
solve  Eqs.  48  and  49  simultaneously  throughout  the  entire  flow 
field.  The  Newton-Raphson  method  (Ralston,  1965)  for  the  inner 
iteration  can  be  expressed  as 


In  Eqs.  48  and  49  the  subscript  notation  has  been  used  to  denote 
the  value  of  the  grid  points  as  follows: 


,j   '     rl    "    ri+l,j  '     r2    ~    ri,j+l   '    r3    =    ri-l,j 


(n+1) 


(n+1) 


(n) 


(n) 


df. 


(n) 


(n) 


Sr 


df 


(n) 


df 


(n) 


(n) 


(n) 


and     r4    =     r.^ 


.(52) 
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where  D  is  the  Jacobian  determinant 
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and  the  superscript  (11)  refers  to  the  iteration  number.  The 
successive  over  relaxation  method  (Forsythe  and  Wasow,  1960)  for 
the  outer  iteration  can  be  expressed  as 
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where  T0  and  ~z~0  are  the  new  values  of  the  variables  r  and  z 
obtained  from  Eq.  52,  co  is  the  over  relaxation  parameter  with  a 
value  between  1  and  2,  and  the  superscript  (m)  refers  to  the  outer 
iteration  number. 

In  the  case  of  Eq.  50  which  involves  only  one  unknown,  a 
simple  form  of  the  Newton-Raphson  iteration  for  nonlinear  equa- 
tions can  be  used  (Kunz.   1957) 
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where  the  superscript  (n)  refers  to  the  iteration  number  and 
f(r0  (n) )  is  the  function  given  by  Eq.  50.  The  value  of  f(r  0(n)  )  will 
be  zero  for  the  correct  value  of  r0  . 

Expressions  for  the  derivatives  in  Eq.  52  can  be  obtained  by 
differentiating  Eqs.  48  and  49  with  respect  to  r0  and  z0 
respectively,  as  given  below 
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Finite  difference  operators 
for  boundary  conditions 

Values  for  the  variable  r  must  be  determined  along  the 
boundaries  of  the  source  surface,  the  ground  surface  and  the  water 
table  as  part  of  thfc  solution  by  satisfying  appropriate  finite 
difference  approximations  to  the  boundary  conditions  already 
described.  Likewise,  values  for  the  variable  z  must  be  established 
along  the  axis  of  symmetry  and  the  cylindrical  outer  boundary. 

Axis  of  symmetry  (7)  -  Q).  The  finite  difference  approximation 
to  the  normal  derivative  condition  dz/di|i  =  0  obtained  by  setting 
the  value  of  z  at  the  nonexistent  grid  point  outside  the  boundary 
equal  to  that  of  its  reflection  within  the  flow  field  and  using  the 
regular  finite  difference  operator  cannot  be  used  at  this  boundary 
since  the  transformation  to  the  <t>ijj  plane  is  not  valid  along  this  line. 
Consequently  the  solution  near  the  axis  is  obtained  on  an  additional 
row  of  grid  points  adjacent  to  the  axis  spaced  one-tenth  of  the 
regular  spacing  from  the  axis  (Jeppson.  1968b).  The  values  of  z  at 
these  added  points  are  denoted  by  zc..  The  derivatives  in  Eq.  27 
with  respect  to  <\i  at  these  added  points  can  be  obtained  by  using 
Taylor's  series  over  an  irregular  grid  network  adjacent  to  the  axis  of 
symmetry  as  shown  in  Fig.  8  in  which  the  dashed  lines  denoted  by  4 
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Fig.  8.      Irregular  grid  network  for  the  special   treatment  of  line 
singularity  at  the  axis  of  symmetry. 
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and  c1  are  obtained  from  the  reflection  of  grid  points  2  and  c  about 
the  axis  of  symmetry.  The  finite  difference  expressions  for  first  and 
second  derivatives  obtained  from  Taylor's  series  (see  Jeppson, 
1968b,  p.  1283)  are 
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To  obtain  a  finite  difference  expression  for  Eq.  27,  the  value  of 
r  at  the  added  grid  point  c  denoted  by  rCj  must  be  known.  It  can 
be  approximated  by  assuming  the  flow  is  uniform  between  the  axis 
of  symmetry  and  the  first  streamline  or  the  first  row  of  regular  grid 
points. 


=     0.  1 
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Eq.  63  is  solved  for  zci  by  the  Newton-Raphson  iterative  method 
(see  Eq.  55).  The  first  derivative  of  f(zCj  )  required  by  the 
Newton-Raphson  equation  is 
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Substituting   Hqs.  60,  61,  and  62  into  Eq.  27,  the  following 
finite  difference  expression  results  for  z  along  the  axis  of  symmetry. 
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The  values  of  z  on  the  axis  of  symmetry  are  obtained  through 
the  use  of  the  finite  difference  expression  for  the  condition  that 
d  z/dO  =  0  in  terms  of  the  values  zn,  zc.  and  zi2.  A  third 
degree  polynomial  approximation  for  z  n    gives 


i,l 


=     1. 01010101   z 


0. 0101010101   z 


i,2 


.(65) 


Source  surface  (7)  -  Q).  A  two  step  operation  is  needed  to 
evaluate  Eq.  35  for  determining  the  boundary  value  of  r.  The  first 
step  is  to  evaluate  the  derivatives  bz/  30,  and  the  second  step  is  to 
perform  numerical  integration  of  Eq.  35.  The  derivatives  have  been 
evaluated  at  each  grid  point  along  the  boundary  by  differentiating 
the  backward  Gregory-Newton  interpolation  formula  (Wylie,  1960) 
truncating  fourth  order  terms.  After  expressing  the  results  in  terms 
of  the  values  of  z  at  the  interior  points 
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Eq.  48,  and  eliminating  the  value  of  r  which  lies  outside  the 
boundary  at  a  nonexistent  grid  point.  The  operator  is  the  same  as 
Eq.  48  with  r4  replacing  r2.  The  two  step  operation  described  in  the 
last  section  can  also  be  applied  to  evaluate  Eq.  36  for  determining 
the  value  of  r  along  this  boundary. 

Vertical  outer  boundary  Q)  -  @.  The  normal  derivative 
condition  Sz/SiJ;  =0  is  replaced  by  second  order  central 
differences  and  the  result  is  combined  with  Eq.  49  to  eliminate  the 
value  of  z  which  lies  outside  the  boundary  at  a  nonexisting  grid 
point.  The  resulting  operator  is  the  same  as  Eq.  49  with  z4  replacing 
z2.  Caution  is  necessary  so  that  a  good  initial  value  of  z  is  supplied 
when  this  is  solved  by  the  Newton-Raphson  method.  Upon 
inspection  of  the  function  f2(r0,  z0)  in  Eq.  48,  the  slope  of  the 
function  at  the  point  where  the  solution  exists  was  found  to  be 
negative  but  changed  to  positive  with  a  small  decrease  in  z.  Before 
the  iteration  begins,  the  slope  of  the  function  with  given  initial  z  is 
examined.  If  the  slope  is  positive,  the  initial  value  of  z  is  increased 
until  the  slope  changes  to  negative.  Subsequently,  the  regular 
iterative  procedure  is  followed  to  obtain  the  solution. 

Water  table\2)  -(J).  To  satisfy  the  normal  condition  3r/  3 $ ,  = 
0  along  the  boundary  a  parabola  is  fit  through  three  consecutive 
points  on  and  adjacent  to  the  boundary  for  each  streamline  as  it 
enters  the  water  table.  The  parabola  is  symmetric  about  the 
boundary.  Two  interior  grid  points  on  the  streamline  adjacent  to  the 
boundary  are  used  to  determine  the  parabola  which  passes  through 
them.  The  radius  at  which  this  parabola  intersects  with  the  water 
table  gives  the  value  of  r  on  the  boundary. 


in  which  NH  is  the  total  number  of  equipotential  lines  specified  for 
the  entire  region,  and  this  boundary  is  specified  as  the  NH  th 
equipotential  line,  and  NH1=NH-1,  NH2=NH-2,  and  NH3=NH-3. 
The  second  step  for  carrying  out  the  numerical  integration  is  based 
on  central  differences  except  at  the  end  grid  points  of  the  boundary. 
The  symmetric  Bessel's  interpolation  formula  is  used  to  evaluate  the 
integration  in  terms  of  two  pairs  of  derivatives  on  each  side  of  the 
grid  point  over  which  the  integration  is  to  be  performed.  Eq.  35 
becomes, 
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Determination  of  the  ratio  R  of  horizontal 
to  vertical  grid  spacing 

The  ratio  of  spacing  between  adjacent  streamlines  to  the 
spacing  between  adjacent  transformed  equipotential  lines  is  R  = 
A  4VA$t.  The  number  of  $t=  constant  lines,  the  number  of  4"  = 
constant  lines  used  in  the  finite  difference  solution,  and  the  total 
head  on  the  source  surface  are  specified  as  input  parameters  in 
defining  a  problem  to  be  solved.  The  value  of  the  ratio  A^/AO,  is 
adjusted  repeatedly  during  the  iterative  process  to  ensure  that  the 
boundary  conditions  of  the  specified  problem  are  satisfied.  To 
adjust  R,  an  expression  for  R  in  terms  of  ijjnand  $s  or  ($,)s 
(where  the  subscripts  s  and  n  may  be  considered  as  curvilinear 
coordinates)  is  needed.  A  relationship  commonly  used  for  the 
general  two-dimensional  case  can  be  expressed  as 


Ai|; 
An 


A<J> 


(68) 


_Ajj_    _    —  An 
A<t>      "         r   As 


Upon  multiplying  both  sides  of  Eq.  66  by  (rNH j+1  +  rNHj),  the 
value  of  rNHj+]  can  be  obtained  when  the  value  of  rNHj  is  known. 
Consequently,  the  integration  of  Eq.  67  begins  at  the  axis  of 
symmetry  and  proceeds  toward  the  outer  edge  of  the  boundary. 
Since  the  values  of  the  derivatives  are  changing  until  the  entire  field 
is  settled,  the  values  of  r  on  the  boundary  must  be  recomputed 
during  each  outer  iteration. 

Ground  surface  ©  -  @.  A  finite  difference  operator  for  r  on 
this  boundary  is  obtained  by  replacing  the  derivative  dr/di|j  =  0 
with  second  order  central  differences,  combining  the  results  with 


where  s  and  n  are  considered  as  orthogonal  curvilinear  coordinates 
commonly  called  natural  coordinates;  s  is  taken  along  any  stream- 
line and  n  is  normal  to  s.  The  spacing  As  is  the  arc  length  along  a  ^ 
=  constant  streamline  between  two  consecutive  $t  =  constant  lines, 
and  An  is  the  arc  length  along  a  <J>t=  constant  line  between 
adjacent  streamlines.  For  the  axisymmetric  case,  consider  a  set  of 
orthogonal  unit  vectors  es  and  en  in  which  the  subscripts  s  and  n 
denote  the  orthogonal  curvilinear  coordinates  as  shown  in  Fig.  9. 
Let  (a\,a2>  and  (a,],a2)  be  the  direction  cosines  of  the  unit  vectors 
es  and  en  respectively.  Then,  since  en  is  a  unit  vector  in  a  normal 
direction  obtained  from  es    by  counterclockwise  rotation  of  90°, 
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/  (W 


«-!•  *2>  * 


Fig.  9.   Orthogonal  unit  vectors. 


and 


chjj 
On 


7^-e        =     <*,  531    +   <  > 

n  1  or  2  Oz 


■  (75) 


Substituting    Eqs.    5,    6,   and   70   into   Eqs.    72    through    75    anc 
rearranging  the  results  gives 


-v  ,  /         oO  a*    \  a* 

OLJJ  /      ,  t  I  t  t 

si  =  Kr  r   di   s7  +  dzs7    =  Kr  r  3r7     •  -(76) 


the  direction  cosines  aj,a2.  ctj,  and  a2  can  all  be  obtained  from  the 
angle  between  the  horizontal  and  the  direction  of  the  streamline  or 

tfj    =     cos   9  ;    d      -    sin  0 

and  (69) 

rtj'=     -   sin  9  ;    d^    =     cos   e 

From  Eq.  69  it  is  obvious  that 

*;  -  -  az 

and  (70) 

*2    =     "l 

Let   er     and   e7     be   the   unit  vectors  in   the  r  and  z  directions, 
respectively.  Then  the  gradient  of  'J1,  can  be  written  as 

od>     __  o$     _^ 

grad  *      =   V$     =    3—    e       +    3—    e  (71) 

B  t  t  or        r  oz        z y      ' 

The  components  of  the  directional  derivatives  of   $  t  and  cp  in  the 
directions  of  es    and  en  are 

0$  _  5$  s*t 

t     =    V*   .  e       ,     d    v-  +  a     v-      (7^) 

Ss  t        s  1   Or  2    dz 

Pi*  a*  0$ 

dQt  —  1         t    .       1         t       (73) 

^     =     V(rt'en     =    M7+a2^ 

|4      =     7+   .T       =    dl|^    +d2|i        (74) 

0s  T         s  1   or  '   oz 


and 


04;  — 


o* 


'  |di  a7  +  d2 


K      r 


An  approximation  of  Eq.  77  is 


An 


A* 
r        As 


Rearranging  Eq.  78  and  solving  for  R  gives 


R    = 


Ai)j 
A* 


—  An 

K      r    — 

r  As 


.(77) 


.(78) 


.(79, 


The  minus  sign  in  Eq.  79  is  due  to  the  fact  that   $ ,  decreases  alonj 
the  positive  s  direction. 

The  arc  lengths  between  the  streamlines  for  each  rectangulai 
mesh  of  the  orthogonal  grid  network  Ant,  and  An2  shown  in  Fig 
1 0  are  approximated  by  straight  lines 


2  2,1/2 

An     (I)    =     [(r.  -  r.    .)     +  (z.  -  z.    .)  .(80) 

1  1,1  +  1  i,j  i.j  +  1  i,J 


and 

An,  (I)    =     [(r. 


-  r  .) 

i+l.j+1  i+l,J 


+  (zi+l,j+l   "  Zi+l,j' 


,2-,  1/2 


(81. 
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'                    m 
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Finally  the  ratio  R  is  obtained  by  averaging  the  R(l)'s  through  two 
selected  consecutive  streamlines 


NH1 


R     =    ^Tl     r^      R(I) 


R     = 


1 


NH1 


NH1 
1  =  1 


K„      (I) 


An   (I)+  An   (I) 
ASl(I)+  As2(I) 


.(87) 


Fig.  10.   Orthogonal  grid  network. 


in  which  the  parenthesized  "I"  denotes  the  Ith  rectangular  mesh 
from  the  lower  boundary  (water  table),  and  the  arc  lengths  between 
the  equipotential  lines  As  j  and  As2  are  approximated  by 


As^I) 


,2       *  *2nl/2 

(ri+l,j-ri,j)    ^-i+lJ-V    ]  .    .(82) 


in  which  NH1=NH-1;  NH  is  the  number  of  equipotential  lines 
specified  as  the  input  data  parameter.  The  ^  =  constant  lines 
between  which  the  ratios  R(I)  are  calculated  are  removed  from  the 
axis  of  symmetry  to  avoid  being  affected  by  the  line  singularity  at 
the  axis  of  symmetry  which  will  be  discussed  in  a  later  section. 

Determination  of  scaled  volume  flux  rate 

In  order  to  solve  Eqs.  48  and  49  simultaneously  as  shown  in 
Eq.  52,  the  scaled  volume  flux  rate  A<\i/K0  must  be  known.  The 
rate  of  flow  between  any  two  consecutive  stream  surfaces  AQ  is 
given  by 


and 


AsJI)  =     [   (r  -  r  )' 

2  l       i  +  l,j  +  l        i,j  +  l 


A$ 
AQ    =     K   —    AA 
As 


+  (z.    .     .        -  z  ) 

i  +  l.j  +  1       i,j+l 


2n   1/2 


.(83) 


in  which  A  is  the  annular  area  perpendicular  to  the  direction  of  flow 
and  surrounded  by  two  consecutive  stream  surfaces.  An  approxima- 
tion to  this  area  using  the  average  radii  and  normal  increments  An  is 


Values  for  Kr  and  r  in  Eq.  79  are  taken  as  the  average  permeability 
parameter  and  average  radial  distance  from  the  axis  of  symmetry  for 
the  four  grid  points  forming  each  mesh.  These  averages  denoted  by 
Kr      and  r     respectively  are  given  by  the  following  equations 


AA    =    -    (r.    .+  r.  +r  .  +  r...    ,,,)(An    +An    ) 

4         i  ,j        i,j+l  l+l  ,j        i+l,j+l  1  2 

(89) 


K  =    -    (K         +K  +K  +K  1 

rm  4         rit-         riJ  +  1         ri+1#j  +Kri+1>j  +  1  J       .(84) 


Substituting  Eq.  89  into  Eq.  88  gives 


in  which  Kr  's  can  be  obtained  from  Eq.  1  7,  and 


T   (r-    •  +  r-    •    ,  +  r-     ,    •  +  r-     ,    •    J 
4        i,j        i,j  +  l        i  +  l.j        i  +  l,j  +  l 


(85) 


AQ  (I)    =    2ir  K       K        (I)  M>  r      (I) 
or  m 

m 


AjijID  +  An^I) 


ASl(I)  +  As2(I) 


(90) 


Thus  the  ratio  A|/A$  tfor  each  rectangular  mesh  is  given  by 


and  the  average  value  of  AQ(I)  is 


R(D     =    KT 


d)rm(I) 


An    (I)+An(I) 
As1  (I)+As2(I) 


(86) 


AQ    = 


NH1 


NH1 
S       AQ(I) 
1=1 


(91) 
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Since  for  three-dimensional  axisymmetric  flow  AQ  =  2tAi4i, 


Subdivided  grid  network 


Aip 
K 


A<J> 


NH1 


S         K        (I)  r      (I) 

NH1  r  mv 

1  =  1  m 


An    (I)+An    (I) 


As    (I)  +  As2(I) 


.(92) 


In  parts  of  the  region  of  the  problem  values  of  Ar  and  Az  are 
greater  than  in  other  regions.  Better  accuracy  can  be  expected  from 
the  finite  difference  solutions  in  regions  directly  underneath  the 
source  circle  where  Ar  and  Az  are  the  smallest.  At  the  region  closer 
to  the  outer  boundary  where  the  moisture  content  is  very  low.  Ar 
and  Az  will  be  large  and  poorer  precision  can  be  expected  in  finite 
difference  solutions.  While  this  outside  region  is  not  of  major 
interest,  the  poor  approximations  there  will  affect  the  accuracy  of 
the  finite  difference  solution  in  other  regions. 


and  from  Eq.  87 


A^j 
K 


A$R      (93) 


However,  if  the  spacing  of  the  grid  network  is  made  too  small, 
unnecessarily  large  amounts  of  computer  time  are  required  for  the 
solution.  The  convergence  rate  of  the  outer  iteration  of  the  system 
also  decreases  with  increasing  number  of  grid  points.  Better  results 
are  obtained  by  subdividing  the  grid  network  over  part  of  the 
region. 


in  which  A<t>  =  H/NH1,  (H  is  the  total  fluid  head  at  the  surface 
inside  the  rainfall  simulator  and  is  specified  as  an  input  parameter). 
Finally  Eq.  93  becomes 


Ai|j 
K 


H 


NH1 


(94) 


and  the  scaled  volume  flux  rate  Av|j/K0   can  be  obtained  from  the 
ratio  R=  Ai|j/A<t>   . 


In  the  present  analysis,  a  finer  grid  network  is  used  for  the 
region  enclosed  by  the  last  two  4*  =  constant  lines  and  the  source 
surface  and  water  table  boundaries.  The  spacing  of  the  subdivided 
grid  network  is  one-quarter  of  the  regular  spacing.  The  values  of  r 
and  z  for  the  grid  points  in  the  subdivided  network  are  related  to 
the  values  at  the  regular  grid  points  as  shown  in  Fig.  1  1.  During  each 
iteration  through  the  regular  grid  network  the  values  of  r  and  z  at 

these  grid  points  represented  by  i  =  1.2,3 NH,  j  =  NS-2  and  i  = 

1,2,3,  ....  NH,  j  =  NS-1  (where  NH  and  NS  represent  the  values  of  i 
and  j  for  the  final  4>t  =  constant  and  ^  =  constant  lines 
respectively)  are  used  to  determine  the  values  of  r  and  z  for  the  grid 


]    =NS2 


Fig.  11.   Subdivided  grid  network  at  one-half  and  one-quarter  of  regular  spacing  along  the  upper 
boundary  of  the  Ot^  plane. 


21 


points  in  the  subdivided  grid  network  denoted  by  i ,  =  1,2,3,  ....  NP, 
j  j  =  1  (where  NP=2NH-1)  with  halt"  of  the  regular  spacing.  A  special 
operator  is  required  to  relate  the  variables  r  and  /  at  the  point  0 
shown  in  Fig.  1  1  with  half  of  the  regular  spacing  to  the  variables  at 
regular  grid  points  1.2.3.  and  4  and  grid  points  5  and  6  on  the 
subdivided  grid  network  with  half  of  the  regular  spacing  (see  Fig. 
11).  Suppose  a  function  f  =  f(  £,,6)  and  a  rectangular  mesh  is 
subdivided  into  halves  of  the  original  spacing  as  shown  in  Fig.  12. 
the  values  of  fat  any  two  points  denoted  by  (i-,.  ^)  and  (  £,+  d.  ^j 
+  m)  in  which  d  =  ';Af  (  and  m  =  ':Ao  are  related  by  the  Taylor's 
series 


f(«t+   d,  i|>+m)    =    £($t,  40+    fd§^  +  m^]     f(°t-  + 


a 


n-l 


(nTT)T     \dWt    +m§^  *(V+,+Rn   •    '(95) 


where  R  n  is  the  remainder  term  and  is  given  by 


and 

I 

o  <  i  <   l 


(96) 


to  indicate  that  the  truncation  of  Rn  from  Eq.  95  has  error  of  order 
(|d|+|m|)n.  To  obtain  the  operator  at  point  0,  Taylor's  series  of  f 
is  written  for  points  1,  2,  3,  and  4  shown  in  Fig.  12  about  the  value 
of  f  at  the  central  point  0, 


dm 


(98a) 


o  \     t '  0  v    T   'o 

32f 
+  dm|  Aj,  ^,  | (98b) 


di|i 


f.     =  f    +d 
3  o 


l;)-m(|)+i 


^-(9 


-  dm 


t 


(98c) 


and 


Equation  96  can  be  rewritten  as 


°  t      n  'O 


Rn     =     0[(|d|    +  |  m|  )n] 


.(97) 


+  dm 


Sdj 


(98d) 


Fig.  1 2.   Subdivision  of  a  rectangular  mesh. 


in  which  third  and  higher  order  terms  are  neglected  and  the  error 
term  is  0{(  |d[+-|m|)3 }  .  The  subscripts  denote  the  points  at  which 
the  function  f  or  its  derivatives  are  taken.  Adding  Eqs.  98a  through 
98d  and  substracting  4  f0  from  both  sides  of  the  equal  sign  of  the 
resulting  equation  gives 


3    f 


n  \      t  '  i 


2 

t  '  < 


2m 


z^i+WV-4^ 


(99) 


Taylor's  series  for  f  5  and  f  6  about  f  „  at  the  central  point  0  can 
also  be  written  to  obtain  a  finite  difference  expression  for  the 
second  partial  derivative  of  f  with  respect  to  ^  t  at  the  central  point. 


a*' 


7  (fs  +  f6 

d 


2f   ) 
o 


(100) 
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By  multiplying  Eqs.  99  and  100  by  Kr2f02  and  {  1  - 
Kr2  f02  (d/m)2},  respectively,  and  combining  the  resulting  equa- 
tions, a  finite  difference  expression  for  the  terms  enclosed  by  the 
first  parenthesis  in  the  first  governing  equation  Eq.  26  at  the  central 
point  is  obtained. 


K2f2    R 

r  °  W2 


So2 


—  2      2 
r     o 


2(fl+f2  +  f3+f4) 


F,(r    ,z    )  =  16(K    )2   r2 
2      o     o  r  o      o 


+  16R    (z    +  z   -2z  ) 
5       6       o 


-(z1+z2  +  z3  +  z4)-(z5  +  z6) 


(k")3   r2   (^f  r,(z   -*    l^'lz+z     -z     -z,)2 


r  o      o     \Pb/  o     o 


12        3        4' 


(P£|n 


+  4R(K    )     f*      r,(z   -S>   ) 
r  o\Pb/  o     o 


n-i    i 


-R(z2+z3-Vz4) 


(f5  +  V] 


+  —  (f .+  f ,  -  2  f    ) 
,2       5       6  o 

d 


(101) 


(z2+z3-zrz4) 


Finite  difference  expressions  for  the  terms  of/d^,     and  S f/oip    at 
the  central  point  can  also  be  obtained  from  Eqs.  98a  through  98d 


+  2R(Kf)o  (z1+z2-z3-z4)(z2+z3-z1-z4)     =     0  (|05) 


*      to 


(102) 


and 


U) 


to  ^i+v-^+v1 


(103) 


Because  the  values  of  r  and  z  at  point  7  of  the  rectangular  mesh 
shown  in  Fig.  1 2  are  not  available  in  this  transition  zone  from  the 
regular  grid  network  to  the  subdivided  grid  network  with  one-half  of 
the  regular  spacing,  the  regular  five-point  operators  Eqs.  48  and  49 
cannot  be  used.  Instead,  the  new  values  of  r  and  z  at  point  0  are 
obtained  by  solving  the  new  seven-point  operators  Eqs.  104  and  105 
simultaneously.  The  Newton-Raphson  iterative  method  is  used  in 
this  inner  iteration  shown  by  Eq.  52,  and  Eq.  54  is  then  applied  for 
outer  iteration.  The  expressions  of  the  derivatives  to  be  used  in  Eq. 
52  are 


Substituting  Eqs.  100,  101,  102,  and  103  into  Eq.  26  with  proper 
substitution  of  r  and  z  for  f,  a  finite  difference  operator  for  Eq.  26 
at  the  central  point  is  obtained.  The  operator  will  be  evaluated  in 
terms  of  the  values  of  r  and  z  at  the  surrounding  grid  points  1.  2,  3, 
4,  5,  and  6  with  half  of  the  regular  spacing.  (See  Fig.  1  2.) 


^1 

T— ~  =  48  (K    ) 

or  r  o 


iww-w 


+  2 


ro\Pb/         o     o 


[ri+r2-r3-r4)(r2  +  r3-rrr4) 


F,  (r    ,z    )  =   16 (K    ) 
loo  r  o 


^(rl  +  r2  +  r3  +  r4)-(r5  +  r6)   |  r^ 


32R2+(K    )  2  (r.+r    -r    -r    )2 

r  o        12       3       4 


r       + 
o 


16R    (r5+r6) 


^Volpff^oV'VvVV'VVVV 


32R2  +  (K    )2   (r.  +r     -  r     -  r    )2 
r  o         1        2        3        4 


16R     (r5  +  r6) 


-«S"r).(f9-«V.),,",(^)<'z+'3 -'.-'«> 


•(r,  +  r    -r    -  r,) 


(104) 


By  making  the  same  substitutions,  a  similar  operator  for  Eq.  27  is 
also  obtained. 


««V.(fg  ^o->/-'(;f)<VVVV 


(106) 


Oz  r  o    \Pi_J  o      o 


t 
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:5+r6-2(rl+r2+r3+r4,|ro3   -(rl  +  r2-r3-r4)2rof 


+  2(pJ)    ^Vr2-Vr4)(r2  +  r3-Vr4)Rro2 

J(F)2(n-D(z  -*/"2  -  2(K-)3(^)\[(Z  -of1]2} 

<j       r  o  o      °  r  o  \Pb/  o       o  J 

-  4R(?]\f1(r?+r    -r  -r    )r     J  (K~ )    (n-l)(z    -O)"1"2 
IP]-,/         K2314oro  oo 


r  o   ^Pb/  o      o  j 


(107) 
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^2 

Z7~ 


16 


|[(zl+z2+z3+z4)-2(z5+z6)] 


'"'l+W*/}    (Vf    ro Il08) 


^   =    -  32  R2  -  ,Kr,o3  r*    K  »<»l+«2-3-4' 


h-ix..-./"'  -3<vopn*vVn"'i2} 


+  R 


(ft)'1 "  [R'V 


z3_zrz4)      '~ 


It) 


(z2+z3-z1-z4)(Kr)o{(r1-l)(zo-$o)^2-(Kr)o(^[(zo-0/-Y} 


-2R(z1+z2-z3-z4)(z2+z3-z1-z4)(Kr)o2(^)r1(zo-$o)T1-1 


(109) 


In  order  to  iterate  over  the  transition  zone  from  the  regular  grid 
network  to  the  subdivided  grid  network,  the  values  of  r  and  z  on  the 
boundary  of  the  subdivided  grid  network  (i.e.,  the  line  denoted  by 
(J!  =  1)  in  Fig.  11)  are  obtained  from  either  the  seven-point 
operators  or  the  five-point  operators  depending  on  whether  the  grid 
point  in  question  is  inside  or  on  the  boundary  of  regular  rectangular 
mesh.  Along  the  boundary  denoted  by  j  j  =  1,  the  values  for  r  and  z 
at   grid   points  which   are   located   at   the,  centers  of  the  regular 

rectangular  meshes  (i.e.,  the  grid  points  denoted  by  (i  l  =  2,4 

NP-1,  Ji  =  1)  in  Fig.  11)  are  obtained  from  the  seven-point 
operators,  Eqs.  104  and  105.  Values  at  the  remaining  grid  points  on 
the  boundary  (j  l  =  1 )  are  obtained  by  solving  the  original  five-point 
operators  with  points  1,  0,  4,  and  7  as  surrounding  points  with  the 
point  5  as  the  point  in  question.  The  order  of  the  iteration  over  this 
boundary  ( j  ,  =  1 )  starts  with  i  j  =  2  and  proceeds  to  i ,  =  MP- 1  ..The 
seven-point  and  five-point  operators  are  used  in  this  last  pass 
alternatively. 

The  values  of  r  and  z  at  the  soil  surface  and  water  table  have  to 
be  approximated  separately  to  satisfy  their  respective  boundary 
conditions.  The  operator  for  the  point  at  the  soil  surface  denoted  by 
(i  i  =  NP.  j  ]  =  1)  is  obtained  by  noting  that  the  value  of  z  is  constant 
at  this  boundary  being  equal  to  the  depth  of  the  water  table.  Since  r 
changes  by  only  a  small  amount  along  the  surface,  the  value  of  r  at 
(i  |  =  MP,  j  [  =  1 )  can  be  taken  as  the  average  of  the  values  of  r  of  the 
two  regular  grid  points  on  the  surface  denoted  by  (i=NH,  j=NS-l) 
and  (i=NH,  j=NS-2)  on  the  regular  grid  network  system.  To  evaluate 
the  value  of  r  at  the  point  denoted  by  (i ,  =  1 ,  j ,  =  1 )  on  the  water 
table,  the  condition  that  the  streamlines  enter  the  water  table  at 
right  angle  (i.e.  d  r/  S  $,  =  0)  must  be  satisfied.  The  value  of  r  can 
be  determined  by  substituting  this  normal  derivative  condition  into 
the  first  governing  equation,  Eq.  48. 

After  computing  values  of  r  and  z  on  the  grid  network  with 
one-half  of  the  regular  spacing  values  of  r  and  z  on  the  grid  network 
with  one-quarter  of  the  regular  spacing  denoted  by  i2  =  1,2,3,  ..., 
MP,  j2  =  1,  where  MP  =  4NH-3  are  computed  in  a  similar  manner. 


After  solving  for  the  values  of  r  and  z  on  the  grid  network  with 
one-quarter  of  the  regular  spacing  by  the  iterative  procedure 
described  in  a  previous  section,  values  of  r  and  /.  at  the  regular  grid 
network  covering  the  subdivided  region  are  equated  to  the  values  of 
r  and  z  at  the  corresponding  grid  points  with  one-quarter  of  the 
regular  spacing. 

Adjustment  of  the  number  of  equipotential  lines 
which  intersect  the  vertical  outer  boundary 

The  number  of  equipotential  lines  which  intersect  the  vertical 
outer  boundary  is  denoted  by  M  4.  If  M4  is  equal  to  or  greater  than 
two,  one  of  these  equipotential  lines  is  assumed  to  pass  through  the 
intersection  of  the  vertical  outer  boundary  and  the  ground  surface. 
The  criterion  used  to  determine  the  value  of  M4  is  based  upon  the 
orthogonality  of  streamline  with  equipotential  line.  Orthogonality 
can  be  checked  at  each  grid  point  by  computing  the  angle  between 
equipotential  line  and  streamline  at  the  point.  If  the  angle  is  90°, 
the  condition  is  satisfied;  otherwise,  M4  needs  to  be  adjusted.  For 
example,  the  angle  ((3)  is  checked  at  the  subdivided  grid  point 
where  the  third  equipotential  line  from  the  water  table  and  the  third 
streamline  from  the  vertical  outer  boundary  intersect.  If  the  angle  is 
90°  ±  1°,  the  orthogonality  condition  is  assumed  to  be  satisfied, 
and  the  value  of  M  4  is  not  changed.  If  the  angle  is  greater  than  91°  , 
M4  is  decreased  by  one,  and  if  the  angle  is  less  than  89°,  M4  is 
decreased  by  one.  The  adjustments  made  for  the  corner  point  and 
the  surrounding  grid  points  are  shown  in  Fig.  1 3. 
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-^%sA^^ 

•     JjtC                                             ' 
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Water    table      | 

(b)        P       >  90  ,    M4  =  M4+l 

Fig.  13.  Adjustments  of  the  number  of  equipotential  lines  inter- 
secting the  vertical  outer  boundary  and  related  grid  points 
in  the  subdivided  grid  network. 
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SOLUTION  AND  DISCUSSION  OF  RESULTS 


Nature  of  solution  and  results 
obtained  from  solution 

Several  solutions  have  been  obtained  for  the  steady-state  flow 
system  resulting  from  moisture  applied  over  a  circular  area  at  the 
ground  surface.  By  varying  the  values  of  the  parameters  in  Eq.  16, 
the  effect  of  each  parameter  on  the  resulting  flow  pattern  is  studied. 
The  solution,  which  is  obtained  through  the  inverse  formulation 
described  earlier,  consists  of  values  of  r  and  z  coordinates  in  the 
radial  and  vertical  directions  respectively  at  each  grid  point  in  the 
finite  difference  solution.  Consequently,  the  solution  output  ob- 
tained can  be  easily  plotted  to  form  a  flownet  to  give  a  general  idea 
about  the  How  pattern.  Flownets  have  been  obtained  from  the 
Gerber  plotter  at  the  University  of  Utah  computer  center  which 
essentially  draws  a  straight  line  between  adjacent  coordinates.  A 
parabolic  interpolation  formula  was  included  in  the  FORTRAN 
plotting  subroutine  to  plot  smooth  curves  in  the  region  near  the 
water  table  where  the  spacings  of  equipotential  lines  are  too  large  to 
obtain  a  good  flownet  by  simply  drawing  a  straight  line  between 
two  adjacent  grid  points. 

In  this  section  a  number  of  flownets  (Figs.  14  through  23)  are 
given  which  represent  the  solutions  obtained  for  varying  soil 
parameters.  All  of  the  solutions  have  specified  that  pc/y  =  1.0  feet 
over  the  source  circle.  These  solutions  can  be  divided  into  three 
groups.  The  solutions  in  the  first  group  (Figs.  14  through  18)  have 
common  pb  and  n  values  (pb/7  =  3.5  feet  and  q  =  2.50)  but 
varying  b  values  ranging  from  1.0  to  3.0.  In  the  second  group,  (Figs. 
18  through  22)  the  solutions  have  common  values  for  b  and  n  (b  = 
1.0  and  n  =  2.50)  and  different  bubbling  pressure  heads  ranging 
from  1.5  feet  and  3.5  feet.  In  the  last  group,  (Figs.  18  and  23)  the 
solutions  have  common  b  and  pb  values  (b  =  1.0  and  pb/7=  3.5 
feet)  and  the  pore  size  distribution  index  q  ranging  from  2.5  to  3.5 
were  specified  to  obtain  solutions. 

The  portion  of  the  flownets  with  dashed  lines  are  the  regions 
where  the  spacings  on  the  $tuj  plane  have  been  reduced  to 
one-quarter  of  the  regular  spacing  to  increase  the  accuracy  of  the 
numerical  approximation  in  these  regions.  Consequently  one-quarter 
as  much  flux  moves  between  dashed  lines  as  between  adjacent  solid 
lines. 

For  each  flownet  obtained,  the  orthogonality  condition  must 
be  satisfied  at  each  grid  point.  For  the  axisymmetric  isotropic 
homogeneous  porous  medium,  the  orthogonality  is  obtained  in  the 
following  manner.  Along  <l>  =  constant  lines 

d  1  d$ 

d$     =    v"    dr  +  ^~  dz     =    ° 
or  dz 

or  (110) 


l£l 


0$/0] 


Along  i\i  =  constant  lines 


ddj     =    v^    dr  +  ^-x    dz 
dr  dz 


=     0 


(111) 


di|;/d: 
d^/di 


After  substituting  Eqs.  5  and  6  into  Eqs.  110  and  111  and 
multiplying  two  equations  together,  the  following  orthogonality 
condition  which  is  independent  of  K  =  K(p)  can  be  obtained  for 
saturated  and  partially  saturated  axisymmetric  isotropic  homo- 
geneous media. 


HI 


Idzj 


(112) 


in  which  the  subscript  i  is  used  to  denote  that  O  remains  constant. 


Coping  with  multiple  roots  of 
finite  difference  operators 

The  finite  difference  operators,  Eqs.  48  and  49,  are  solved 
simultaneously  by  the  Newton-Raphson  method  at  each  grid  point. 
At  most  of  the  interior  grid  points,  rough  initial  estimates  for  r  and 
z  are  adequate  to  insure  that  the  Newton-Raphson  method  will 
select  the  physically  correct  roots  from  the  operators.  At  some  grid 
points,  particularly  those  close  to  the  vertical  outer  boundary  and 
the  water  table,  this  is  not  the  case.  The  multiple  possible  roots 
resulting  from  simultaneously  solving  Eqs.  48  and  49  are  extremely 
difficult  to  analyze.  For  instance,  in  obtaining  solution  for  q  = 
2.50,  Pb/7  =  1.50 ',  and  b  =  1.0,  it  was  discovered  that  a  slight 
variation  of  the  initial  values  of  z  at  the  grid  points  in  the  vicinity  of 
the  vertical  outer  boundary  and  the  water  table  lead  to  an 
alternative  solution.  Tables  1  and  2  demonstrate  the  sensitivity  of 
the  solution  to  the  initialization.  In  Tables  1  and  2,  f,  is  the 
operator  Eq.  48,  f2  is  the  operator  Eq.  49,  n  is  the  number  of  the 
Newton-Raphson  iteration,  and  r0  and  z0  are  the  initial  values  of  r 
and  z  used  in  obtaining  the  solution.  In  this  case,  starting  with  r0  = 
25.906  and  z0  =  5.0  leads  to  a  final  r  =  25.906  and  z  =  5.295 
whereas  starting  with  r0  =  25.906  and  z0  =  4.5  results  in  the 
solution  r  =  25.894  and  z  =  4.1  15.  When  two  roots  are  relatively 
close  to  each  other  as  shown  in  the  example,  it  is  difficult  to  pick  up 
the  correct  solution  as  required  by  the  physical  system.  Only  after 
the  iterative  solution  has  been  obtained,  and  the  orthogonality 
condition  is  not  satisfied,  is  it  apparent  that  an  incorrect  root  was 
selected  at  some  stage  during  the  solution  process.  The  best 
approach  seems  to  be  to  modify  the  initialization  in  another 
attempt  at  a  solution  in  hopes  that  the  correct  roots  to  the 
functional  operators  will  be  supplied  by  the  Newton-Raphson 
iteration  throughout  the  entire  solution  process. 
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Fig.  20.  Flownet  for  axisymmetric  infiltration  from  source  circle  to  water  table  with  pb/7  =  2.50  ',  b  =  1 .00,  t]  =  2.50,  and  pc/7  =  1 .0' 
over  the  source  circle. 


Fig.  21.  Flownet  for  axisymmetric  infiltration  from  source  circle  to  water  table  with  pb/7  =  2.00  ',  b  =  1.00,  ^  =  2.50,  and  pc/7  -  1.0' 
over  the  source  circle. 


Fig.  22.  Flownet  for  axisymmetric  infiltration  from  source  circle  to  water  table  with  pb/7  =  1.50  ',  b  =  1.00,  *1  =  2.50,  and  p  /7  =  1.0' 
over  the  source  circle. 
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Fig.  23.  Flownet  for  axisymmetric  infiltration  from  source  circle  to  water  table  with  pb/7  =  3.50  ',  b  =  1 .00,  *)  =  3.00,  and  p  /y  =  1 .0' 
over  the  source  circle. 
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Table  1 .   Newton-Raphson  solution  for  rQ  =  25.906    and  z0  =  5.000' . 


n 

r 

z 

df 
l 

Sf2 

bi2 
Sz 

f] 

h 

0 

25.906 

5.000 

1 

25.922 

5.431 

-404.043 

63.019 

-14.505 

-125.209 

-20.630 

54.159 

2 

25.908 

5.307 

-373.712 

77.449 

-11.121 

-295.387 

4.330 

-36.586 

3 

25.906 

5.295 

-380.721 

75.255 

-1  1.944 

-246.546 

0.166 

-  3.011 

4 

25.906 

5.295 

-381.435 

74.907 

-12.030 

-241.772 

0.003 

-  0.029 

5 

25.906 

5.295 

-381.442 

74.903 

-12.030 

-241.725 

0.000 

0.000 

Table  2.   Newton-Raphson  solution  for  rQ  =  25.906    and  zQ  =  4.50' . 


n 

r 

z 

dfl 

*1 

dz 

Sf2 

az 

f! 

h 

0  . 

25.906 

4.500 

1 

25.844 

3.680 

-463.1  10 

-   10.824 

-20.865 

81.981 

-  37.697 

65.909 

2 

25.823 

3.995 

-674.946 

-673.458 

-42.694 

539.558 

198.263 

-170.731 

3 

25.886 

4.103 

-563.347 

-287.334 

-31.833 

333.134 

66.455 

-  33.944 

4 

25.894 

4.115 

-540.024 

-188.164 

-28.999 

274.494 

6.316 

-     3.108 

5 

25.894 

4.115 

-537.621 

-178.362 

-28.702 

268.159 

0.069 

-     0.037 

6 

25.894 

4.115 

-537.591 

-178.248 

-28.698 

268.082 

0.000 

0.000 

Determination  of  final  radius 
of  boundary  value  problem 

The  total  number  of  the  streamlines  of  the  t,^  plane,  NS,  is 
specified  as  an  input  parameter  in  defining  the  problem.  The  value 
NS  =  19  has  been  used  to  obtain  all  solutions  given  herein.  The  NS  = 
19  stream  surface  consists  of  the  ground  surface  outward  from  the 
edge  of  the  source  circle  and  the  vertical  outer  cylindrical  boundary 
of  the  flow  system.  The  second  stream  surface  from  the  vertical 
outer  boundary  in  the  regular  grid  network  system  is  the  stream 
surface  enclosing  95  percent  of  the  total  flux. 


The  maximum  radius  of  the  95  percent  total  flux  stream 
surface  of  the  first  group  of  solutions  has  been  plotted  as  a  function 
of  the  value  of  the  parameter  b  in  Fig.  24.  The  plot  shows  that  the 
radius  increases  linearly  with  b  (pb/7  and  t\  are  held  constant  at 
3.50  feet  and  2.50  respectively).  This  general  relationship  might 
have  been  anticipated  from  Fig.  4,  since  for  a  constant  capillary 
pressure  head,  the  height  of  the  plateau  or  the  maximum  value  of 
K~  increases  as  the  value  of  b  decreases.  That  is,  if  the  two 
parameters  p  /y  and  t)  on  the  source  circle  are  held  constant,  the 
soil  will  become  more  saturated  and  the  flow  system  will  become 
nearly  one-dimensional  with  decreasing  value  of  b. 


4.0 


3.0  - 


2.0- 


3.0 


r  /  R, 


Fig.  24.  The  radius  of  95  percent  flux  stream  surface  at  the  water  table 
as  a  function  of  b  values  with  pb/7  =  3.50' ,  and  r\  =    2.50. 
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Fig.  25.  The   radius   of  95  percent  flux  stream  surface  at   the  water  table  as  a 
function  of  pb/y  values  with  b  =  1 .0.  and  1  =  2.50. 


The  relationship  is  also  linear  between  the  radius  of  the  95 
percent  flux  stream  surface  at  the  water  table  and  the  bubbling 
pressure  head  as  shown  in  Fig.  25.  The  maximum  95  percent  radius 
increases  as  the  bubbling  pressure  head  increases.  Since  the  bubbling 
pressure  increases  with  decreasing  average  pore  size  of  the  soil,  the 
soil  throughout  the  region  becomes  more  saturated  as  pb/7 
decreases.  Thus,  the  flow  system  of  the  problem  becomes  more 
nearly  one-dimensional  as  the  bubbling  pressure  decreases  and  the 
remaining  soil  parameters  are  held  constant.  That  is  the  radius  of  the 
95  percent  flux  stream  surface  at  the  water  table  decreases  as  the 
bubbling  pressure  decreases  as  shown  in  Fig.  25. 

The  region  within  the  95  percent  flux  stream  surface  obtained 
from  the  solutions  for  varying  r\  values  with  pc/7  =  1.0'  specified 
on  the  source  surface  has  the  relative  permeability  close  to  unity. 
Thus,  within  this  region,  the  saturation  increases  as  the  values  of  r\ 
increases  (Fig.  5).  and  the  radius  of  95  percent  flux  stream  surface 
at  the  water  table  decreases  a  relatively  small  amount  as  the  value  of 
r|  increases  (Fig.  26).  The  indicated  tendency  may  not  be  true  when 
a  high  constant  capillary  pressure  is  specified  on  the  source  surface. 
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Determination  of  seepage  velocity 

Expressions  for  computing  the  seepage  velocity.  V,  at  any 
point  can  be  obtained  by  defining  the  Jacobian  for  the  transforma- 
tion of  variables.  Eq.  30  can  be  rewritten  as 


J    = 


Sr      dz 
Suj     Sdj 


(113) 


Substituting  Eqs.  5  and  6  into  Eq.  1  13  gives 
V 


J     = 


V 


K  K 


r  V  r  V 


r  V 
K 


(114) 


The  derivatives  of  <J>  and  ^  with  respect  to  r  and  z  are  related  to  the 
derivatives  of  r  and  z  with  respect  to  $  and  4*  by  J.  Substituting  Eq. 
1 14  into  these  derivatives,  the  following  expressions  are  obtained. 


K  cos  8 

V 


(115) 


Fig.  26.  The  radius  of  95  percent  flux  stream  surface  at  the  water 
table  as  a  function  of  n.  values  with  pbAy  =  3.50' ,  and  b  = 
1.0. 
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in  which  6  is  the  angle  between  the  direction  of  flow  and  the 
horizontal  r-axis.  By  eliminating  the  angle,  any  of  the  following  four 
equations  results  to  obtain  the  velocity. 


By  numerically  approximating  the  derivatives  in  Eqs.  1  19,  120,  121, 
or  122,  the  seepage  velocity  at  any  grid  point  can  be  obtained.  An 
example  of  the  seepage  velocity  distribution  is  shown  in  Fig.  27(a). 
The  time  required  by  the  water  to  travel  from  the  source  surface  at 
a  distance  r/R0  from  the  axis  of  symmetry  to  the  water  table  for 
Pb/"y  =  3.50',  r|  =  2.50,  b  =  1.0  and  D  =  16'  is  shown  in  Figure 
27(b). 

Hydrostatic  pressure  distribution 
and  pressure  gradient 

The  hydrostatic  pressure  or  pore  pressure  at  any  point  can  be 
obtained  from  the  solution  output  of  z(  <t>t,  *\>)  with  the  equation 


V    = 


(!M 


Bi)j| 


(119) 


pgf* 


7  [(I-  1)A$  -  z] 


(123) 


V   =    - 


arr    /a 


(120) 


Since  a  constant  capillary  pressure  head  is  specified  on  the 
surface  of  the  circular  rainfall  simulator  and  there  is  no  other  source 
of  moisture  into  the  system,  the  entire  region  will  be  partially 
saturated,  and  the  pressure  at  each  point  will  be  negative.  By 
definition  this  negative  pressure  is  a  positive  capillary  pressure,  pc. 
Thus,  Eq.  123  can  be  rewritten  as 


2  la*)  +  r    W 


V    = 


1 


2  W   +  r     a+J 


K 


(121) 


(122) 


=    7[z-(I  -  1)A$] 


(124) 


Fig.  28  gives  the  capillary  pressure  distribution  obtained  from  a 
solution  for  pb/7  =  3.50 ',  b  =  1.00,  and  r|  =  2.50  with  a  unit 
capillary  pressure  head  on  the  surface  of  the  source  circle.  The 
capillary  pressure  distribution  of  the  solution  is  similar  to  the 
distribution  of  permeability.  Analysis  of  the  distribution  of  perme- 
ability for  different  solutions  will  be  made  in  a  subsequent  section 
of  this  study. 
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PC/,  =1-0 

llllll.M    I   111 


Fig.  27(a).    The  seepage  velocity  distribution  for  pb/7  =  3.501 ,  b  =  1 .0.  and  n  =  2.50. 
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Fig.  27(b).  Time  (T  in  hour)  needed  by  water  to  travel  from  the  source  surface  at  a  distance  r/R0 
from  the  axis  of  symmetry  to  the  water  table.  <pb/7  =  3.5'  ,  r\  =  2.50,  b  =  1 .0,  n  = 
porosity,  and  Ks  =  permeability  at  saturation  in  fps.) 


Fig.  28.   Capillary  pressure  distribution  for  pb/7  =  3.50 ',  b  =  1.0.  and  1  =  2.50. 
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The  pressure  gradients  can  be  obtained  by  dividing  the  pressure 
difference  between  adjacent  grid  points  by  the  grid  spacing. 


Distribution  of  the  relative  infiltration 
rate  (lr  )  over  source  circle 

Infiltration  will  be  expressed  in  dimensionless  terms  by 
dividing  by  the  permeability  K0.  This  dimensionless  quantity  lr  = 
I/K0  (where  I  is  the  infiltration  rate  per  unit  area)  will  be  referred  to 
as  the  relative  infiltration  rate.  Distribution  of  relative  infiltration 
rate  over  the  source  circle  for  several  values  of  bubbling  pressure  and 
r\  =  2.50,  b  =  1.0,  and  pc/7  =  1.0'  over  the  source  circle  are  shown 
in  Fig.  29.  These  distributions  were  obtained  from  the  mean  relative 
infiltration  rate  between  each  pair  of  consecutive  stream  surfaces 
(the  infiltration  flux  between  stream  surfaces,  2iTAkJj/K0  divided  by 
its  area).  The  infiltration  rate  is  nearly  constant  except  near  the  edge 
of  the  source  circle  where  rapid  increase  of  I  r  occurs. 

The  distributions  of  relative  infiltration  rate  over  the  source 
circle  for  different  b  values  with  pb/7  =  3.50',  r\  =  2.50,  and  pc/7 
=  1.0'  are  given  in  Fig.  30.  The  influence  of  edge  effect  becomes 
less  as  the  value  of  b  decreases. 


The  distribution  of  the  relative  infiltration  rate  over  the  source 
circle  due  to  different  pore-size  distribution  indexes  is  shown  in  Fig. 
31.  From  the  range  of  n. -values  for  which  solutions  have  been 
obtained  it  appears  that  the  distribution  of  the  relative  infiltration 
rate  over  the  source  circle  is  not  very  sensitive  to  changes  in  n.  This 
lack  of  sensitivity  might  be  expected  because  the  functional 
relationship  between  the  permeability  and  the  capillary  pressure 
head  is  insensitive  to  changes  in  n  values  at  pc/7  =  1.0'  as  shown  in 
Fig.  5. 

Distribution  of  permeability  parameter  (Kr)  on 
the  ground  surface  beyond  the  source  circle 

In  partially  saturated  moisture  movement  it  is  useful  to  have 
information  regarding  saturation  or  effective  permeability  through- 
out the  flow  region.  In  this  and  subsequent  sections  the  relation- 
ships of  Kr  with  the  soil  parameters,  as  obtained  from  analyses  of 
the  solution  results,  are  presented.  Since  the  value  of  Kr  at  any  grid 
point  of  the  problem  can  be  obtained  directly  from  the  computer 
solution  results,  and  it  would  require  considerable  computer 
execution  time  to  solve  Eq.  44  for  Se  from  a  known  Kr  or  Kr  at 
each  grid  point,  the  distributions  of  Se  on  the  ground  surface 
beyond  the  source  circle  as  well  as  any  other  surface  or  plane  of 


Fig.  29.  Distribution  of  relative  infiltration  rate,  I  r,  over  source  circle  as  a  function 
of  distance  from  the  center  for  different  Pt,/Y  values  with  r|  =  2.50,  b  = 
1.0,  K0  =  0.001  fps,  andpc/7  =  1.0'. 
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interest  is  not  given  directly  in  this  report.  Instead,  the  distributions 
of  Kr  are  given  and  discussed  in  this  and  subsequent  sections.  The 
distributions  of  Sc  can  be  obtained  from  the  Kr  vs.  Se  relationships, 
Fig.  6  or  7,  as  discussed  earlier.  One  need  only  enter  these  figures 
with  Kr  and  the  appropriate  soil  parameters  b,  r\  and  pb/7  and 
obtain  from  the  figures  the  corresponding  value  of  Se. 

The  distribution  of  permeability  on  the  ground  surface  beyond 
the  circle  of  application  is  quite  different  from  that  within  the 
source  circle.  The  effects  of  each  of  the  parameters  in  Eq.  1  on  the 
distributions  of  Kr  beyond  the  source  circle  are  given  in  Figs.  32, 
33,  and  34. 

Fig.  32  shows  the  series  of  lines  result  from  solutions  based  on 
different  values  of  pb/7    but  with  other  parameters  in  Eq.  1  held 


constant  (  r\  =  2.50,  b  =  1.0).  Since  in  all  of  these  solutions  the 
capillary  pressure,  specified  on  the  source  circle,  was  less  than  the 
bubbling  pressure  each  of  the  curves  begins  at  Kr  =  l/[(  -y/pb)2'5  + 
1]  and  decreases  with  the  distance  beyond  the  source  circle.  For 
smaller  values  of  pb/7  the  permeability  parameter  decreases  more 
rapidly  at  small  radial  distances  than  for  larger  values  of  pb/7-  At 
larger  radial  distances  the  percentage  decrease  in  Kr  is  less  for  small 
values  of  pb/7- 

The  variation  of  the  permeability  parameter  (K, )  with  radial 
distance  from  the  solutions  with  different  values  of  b  specified  but 
with  r\  and  pb/~y  held  constant  have  been  plotted  in  Fig.  33.  In 
contrast  with  the  reduction  in  the  permeability  parameter  (Kr)  with 
decreasing  values  of  pbAy  shown  in  Fig.  32,  Fig.  33  indicates  that 
the   permeability    parameter  (Kr )   decreases   as   the   value  of  the 
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Fig.  32.   Distribution  of  permeability  parameter,  Kr ,  on  the  surface  as  a  function  of  distance  from 
the  edge  of  the  source  circle  for  different  pb/7  values,  with  T  =  2.50,  and  b  =  1 .0. 
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Fig.  33.    Distribution  of  permeability  parameter.  Kt .  on  the  surface  as  a  function  of  distance  from 
the  edge  of  the  source  circle  for  different  b  values,  with  1  =  2.50,  and  pb/7  =  3.50 ' . 


parameter  b  is  increased.  The  percentage  reduction  in  Kr  with  large 
values  of  b  is  greater  at  smaller  radial  distances  that  at  larger  radial 
distances. 


The  manner  in  which  the  permeability  parameter  (Kr  )  varies 
with  radial  distance  beyond  the  source  circle  for  those  solutions 
with  different  values  of  r\  specified  is  shown  in  Fig.  34.  At  small 
radial  distances  (r/R  0  <  1 . 1 )  the  permeability  parameter  ( K  r )  shows 
a  slight  increase  with  larger  specified  values  of  n,.  At  large  radial 
distances,  however,  the  permeability  parameter  (Kr)  decreases 
considerably  with  increasing  values  of  n, .  The  latter  relationship 
occurs  in  regions  where  the  capillary  pressure  is  large,  and  the 
former  where  the  capillary  pressure  is  small.  The  variation  of  Kr 
with  i-|  is  closely  related  to  the  relationship  shown  in  Fig.  5  in  which 


the  separate  curves  all  passed  through  the  point  pc  =  pb  and  Kr  = 
l/(l+b). 

Distributions  of  permeability  parameter 
(Kr  )  along  the  axis  of  symmetry 

The  distributions  of  permeability  along  the  axis  of  symmetry 
resulting  from  the  solution  obtained  for  different  values  of  pb/7 
with  r|  =  2.50  and  b  =  1.0  are  shown  in  Fig.  35.  Since  all  of  these 
solutions  specified  pc/7  =  1.0'  on  the  surface  within  the  source 
circle,  the  permeability  (Kr)  at  the  surface  is  given  by_  Kr  = 
l/[(7/pb'n  +bl-  At  the  water  table  pc  =  0  and  therefore  Kr  =  1/b 
for  any  value  of  pb.  The  results  plotted  in  Fig.  35  show  that  the 
variation  of  the  permeability  is  large  just  above  the  water  table  for 
small  values  of  pb/-y  . 
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Fig.  34.   Distribution  of  permeability  parameter,  K, ,  on  the  surface  as  a  function  of  distance  from 
the  edge  of  the  source  circle  for  different  T  values,  with  pb/7  =  3.50 ' ,  and  b  =  1 .0. 


The  distribution  of  the  permeability  along  the  axis  of 
symmetry  for  different  values  of  b  with  r|  =  2.50,  and  pb/7  =  3. 50" 
is  given  in  Fig.  36.  At  the  water  table  the  permeability  changes  with 
b  since  K]  =  1/b,  but  this  is  only  because  of  the  definition  of  Kr . 
The  permeability  at  the  surface  is  given  by  K^  =  l/KT/Pb)11  +  bl  ■ 
The  quantity  (7/pb)  (where  pb  and  r\  are  given  constants)  in  the 
denominator  becomes  less  significant  as  b  increases.  Thus,  the 
permeability  becomes  more  nearly  constant  along  the  axis  for  larger 
values  of  b. 

The  distributions  of  the  permeability  along  the  axis  of 
symmetry  as  given  by  solutions  for  different  values  of  r)  with  pb/7 
and  b  constant  are  given  in  Fig.  37.  The  relative  permeability  at  the 
surface  is  given  by  Kr  =  l/l^/Pb^  +  b] .  As  r|  increases,  Kr 
increases  approaching  a  constant  value  equal  to  1/b.  An  increase  of 
the  pore-size  distribution  index  indicates  that  the  porous  media  is 


becoming  coarser  and  the  resulting  flow  system  becoming  more 
nearly  one-dimensional. 


Distribution  of  equi-permeability  curves 

The  variation  of  the  effective  saturation,  Se,  within  the  flow 
region  under_  study  is  much  less  sensitive  than  the  permeability 
parameter,  Kr.  The  permeability  parameter  (Kr )  varies  in  a  wider 
range  from  0.1  to  1.0  while  Se  varies  about  from  0.6  to  1.0.  In 
order  to  show  a  clearer  picture  of  the  flow  system,  the  distribution 
of  equi-permeability  lines  (which  are  lines  of  constant,  Kr)  is 
presented  in  this  section.  Actually  these  distributions  may  also  be 
converted  to  distribution  of  equieffective  saturation  lines  from  the 
solutions  of  Eq.  45,  as  shown  in  Figs.  6  and  7. 
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Fig.  35.   Distribution  of  permeability  parameter,  Kr ,  along  the  axis  Fig.  36.   Distribution  of  permeability  parameter,  Kt ,  along  the  axis 

of  symmetry  for  different  pb/7  values,  with  1  =  2.50,  b  =  of  symmetry  for  different  b  values,  with  pb/7  =  3.50'  ,  r\ 

1.0,  and  pc/7  =  1.0  '  over  the  source  circle.  =  2.50,  and  pc/7  =  1.0'  over  the  source  circle. 


The  distribution  of  equi-permeability  lines  throughout  the 
seepage  region  can  be  obtained  by  substituting  the  solution  of  z  into 
Eq.  123  to  compute  the  pressure,  which  in  turn  is  used  in  Eq.  1  or 
Eq.  16  to  compute  Kr.  The  equi-permeability  lines  are  obtained  by 
interpolations  between  the  values  for  Kr  on  the  Figs.  14  through 
23.  The  distribution  of  equi-permeability  lines  for  the  cases  studied 
are  shown  in  Figs.  38  through  48.  How  these  distributions  are 
related  to  the  soil  parameters  in  Eq.  1  is  presented  in  the  following 
section. 

The  influence  of  soil  parameters  on 
the  distribution  of  permeability 

Changes  in  the  Kr  =  0.30  permeability  curve  with  different 
values  of  pb/7  and  1  =  2.50,  and  b  =  1.0  are  given  in  Fig.  49.  The 
curves  are  shifted  outward  from  the  axis  of  symmetry  as  the  value 
of  pb/7  increases,  and  the  portion  of  the  curves  closer  to  the  water 
table,  where  the  capillary  tension  is  high,  has  been  shifted  more, 
since  the  variation  of  the  permeability  parameter,  Kr,  is  greater  at 
larger  bubbling  pressure  heads  as  shown  in  Fig.  3. 

The  variation  of  Kr  =  0.30  curve  for  different  values  of  b  with 
r\  =  2.50,  and  pb/"/  =  3.50'  is  given  in  Fig.  50.  The  curves  are 
shifted  toward  the  axis  of  symmetry  as  the  value  of  b  increases.  The 
shifting  is  more  significant  in  the  area  of  low  capillary  tension 
underneath  the  source  surface  where  the  variation  of  Kr  is  more 
significant  for  varying  value  of  b  as  shown  in  Fig.  4. 


The  equi-permeability  curves  for  Kr  =  0.20,  0.30,  0.50,  and 
0.80  for  different  values  of  r\  are  given  in  Fig.  51.  The  permeability 
(Kr )  varies  only  slightly  on  the  source  surface  where  pc  ly  =  1.0  is 
specified  (Fig.  5),  and  it  appears  that  the  Kr  =  0.80  curve  or  any 
other  equi-permeability  curve  with  Kr  greater  than  1/(1  +  b)  =  0.50 
is  shifted  outward  as  the  value  of  r\  increases  and  the  equi- 
permeability  curves  with  Kr  equal  or  less  than  0.50  are  shifted 
inward  as  the  value  of  r\  increases  and  the  shifting  is  more 
significant  in  areas  of  higher  capillary  tension  such  as  the  vertical 
outer  boundary. 

These  same  variations  with  n,  are  revealed  in  F_ig.  4  which 
shows  the  variation  of  the  permeability  parameter,  Kr,  with  the 
capillary  pressure  for  different  pore-size  distribution  indexes. 

Effects  of  soil  parameters  on 
the  "spreading  effect" 

The  lateral  movement  of  flow  in  soil  from  an  infiltrometer  or  a 
source  circle  will  be  affected  by  the  size  of  the  infiltrometer  or  the 
source  circle,  the  surface  head,  and  the  soil  properties  defined  by 
Eq.  1.  Marshall  and  Stirk  (1950)  studied  the  effect  of  lateral 
movement  of  water  in  soil  from  flooded  plots  and  found  that  the 
relative  lateral  movement  is  greatest  in  the  smallest  plot.  They 
introduced  a  correction  factor  to  compensate  for  lateral  movement 
and  reduce  the  influence  of  plot  size.  Schiff  (1953)  studied  the 
effect   of  surface   head   on    final    infiltration    rates   based   on   the 
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Fig.  37.  Distribution  of  permeability  parameter,  Kf ,  along  the  axis 
of  symmetry  for  different  n  values,  with  pb/7  =  3.50'  ,  b 
=  1 .0,  and  pc/7  =  1 .0*  over  the  source  circle. 


performance  of  ring  infiltrometers  and  ponds  and  concluded  that 
final  infiltration  rates  increased  directly  proportional  to  the  head  in 
the  infiltrometer.  A  purpose  of  this  study  is  to  determine  the  effect 
of  each  soil  parameter  in  Eq.  1  on  the  lateral  movement  of  moisture 
in  soils  from  moisture  supplied  over  a  circular  source  area  with 
radius  R  0  =  10.0 '  toward  a  water  table  at  a  depth  of  D  =  1 6.0 ' .  The 
magnitude  of  the  "spreading  effect"  can  be  characterized  by 
dividing  the  three-dimensional  relative  infiltration  rate,  13,  by  the 


one-dimensional    relative    infiltration    rate, 
corresponding  one-dimensional  solution. 


I  j,    obtained    from    a 


The  relative  infiltration  rate  for  one-dimensional  steady  down- 
ward flow  to  a  water  table  can  be  obtained  from  the  capillary 
pressure  distribution  equation  during  steady  downward  flow  toward 
a  water  table.  Arbhabhirama  and  Kridakom  (1968)  used  an 
equation  similar  to  Eq.  1  proposed  by  King  (1965)  except  b=  1.0  to 
obtain  integral  form  of  the  Scott-Corey  (1961)  equation, 


z     =      yPc    7(dpc/Yl/[l-(q/Kr)] (125) 

where  z  is  the  elevation  above  the  water  table,  q  is  the  volume  flow 
per  unit  area  in  direction  z,  and  obtained  a  simple  equation  to 
represent  the  capillary  pressure  distribution  during  steady  one- 
dimensional  downward  flow.  If  Eq.  1  is  substituted  into  Eq.  125 
and  the  integration  is  performed  by  using  the  convergent  series 
suggested  by  Arbhabhirama  and  Kridakorn  (1968),  the  following 
expression  for  the  capillary  pressure  distribution  during  steady 
downward  flow  can  be  obtained  (see  Appendix  A). 


1  -bq 


1  - 


n  +  1 


In 


q  p 

^0*0 
1  -bq. 


(126) 


Fig.  38.   Distribution  of  equi-permeability  curves  for  n  =  3.50,  Pb/7  =  3.50' ,  and  b  =  1 .0. 
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Fig.  51 .  The  Kf  =  0.20,  0.30,  0.50,  and  0.80  equi-permeability  curves  obtained  from  solutions  for 
different  r\  values,  with  pb/7  =  3.50'  and  b  =  1.0. 


where  n  >  2,  z„  =  z/(pb/y),  p0  =  pc/pb,  and  qQ  =  q/K0.  For  the 
region  where  p0  <  1  and  the  values  of  p^  are  small,  Eq.  126  can 
be  approximated  as  the  straight  line. 


h  Pb 

—  =  0.  66096  -  3.  54250  —  +  0.  91072  log  b 

+  0.  46874  log  r)    (i28) 


1  -bqr 


(127) 


If  the  soil  parameters,  the  elevation  above  the  water  table,  and 
capillary  pressure  defined  over  the  surface  are  known,  the  relative 
infiltration  rate  per  unit  area  q„  can  be  obtained  from  Eq.  126.  and 
the  total  intake  of  water  I  j  from  the  one-dimensional  case  which 
corresponds  to  the  flux  1  3  in  the  axisymmetric  case  is  obtained  by 
multiplying  q0  by  the  area  of  the  source  circle.  The  relative 
infiltration  rate  for  the  axisymmetric  case,  I3,  has  been  obtained 
from  the  computer  solution  results. 

Values  for  I3/I1  are  plotted  in  Fig.  52  as  a  function  of  the 
scaled  bubbling  pressure  head,  pb/(7  D),  with  r\  =  2.50,  b  =  1 .0,  and 
D  =  16.0'  (the  depth  of  the  water  table  below  the  source  circle). 
This  plot  shows  that  I3/I1  is  directly  proportional  to  the  scaled 
bubbling  pressure  head.  1 3/I  1  increases  as  the  dimensionless  bubbl- 
ing pressure  head  pb/(-yD)  increases.  Fig.  53  shows  1 3/I  x  plotted 
against  b  with  pb/"y  =  3.50',  and  r\  =  2.50.  This  relationship  is  a 
straight  line  on  the  semi-log  plot,  with  I3/I1  increasing  as  b 
increases.  Fig.  54  contains  a  plot  of  I3/I1  against  the  pore-size 
distribution  index  rj  .  A  straight  line  on  the  semi-log  paper  also 
defines  this  relation. 

In  order  to  define  the  total  relationship  between  the  soil 
parameters  and  the  lateral  flow  of  moisture,  a  normal  regression 
analysis  which  minimizes  the  sum  of  squares  of  deviations  in  the 
normal  direction  (Jeppson  and  Huber,  1969)  was  performed  to 
arrive  at  a  simple  expression  relates  I3/I1  to  the  soil  parameters  in 
Eq.  1.  The  resulting  regression  equation  is 


0.25 


Fig.  52.   Scaled  relative  infiltration  rate  (I3  /1 1 )  as  a  function  of 

scaled  bubbling  head  (pb/(7D)),  with  r)  =  2.50,  and  b  =  1.0. 
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in  which  the  correlation  coefficient  is  0.9995.  A  nomogram  was 
constructed  from  Eq.  128  and  provides  a  relationship  for  I  ,/l  ,  as  a 
function  of  rj,  b,  and  pb/(  ,  D).   this  nomogram  is  given  as  fig.  55. 

Values  of  1 3/I  t  can  be  considered  a  correction  factor  in 
obtaining  <one-dimensional  infiltration  rates  from  circular  infiltro- 
meter  measurements.  Let 


c<:  -> 


then 


(129) 


2    2 


30  40 


(130) 


For  any  known  soil  (i.e.  the  soil  parameters  in  Eq.  1  are  known)  C( 
can  be  obtained  from  the  nomogram  Fig.  55,  and  one-dimensional 
relative  infiltration  rate,  I  j,  can  be  obtained  from  Eq.  130. 


I   8 


^    1.6    - 


4.0         5.0 


Fig.  53.   Scaled  relative  infiltration  rate  < 1 3/I  ( )  as  a  function 
of  b,  with  pb/7   =  3.50',  and  n   =  2.50. 


Fig.  54.  Scaled  relative  infiltration  rate  ( 1 3/I  | )  as  a  function  of  ti, 
with  pb/7  =3.50',  and  b=  1.0. 
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Fig.  55.   Nomogram  which  relates  the  ratio  of  the  axisymmetric  infiltration  rate  divided  by  the 
one-dimensional  rate  to  the  soil  parameters. 
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SUMMARY  AND  CONCLUSIONS 


Solutions  for  the  problem  of  steady-state  partially  saturated 
axisymmetric  infiltration  resulting  from  moisture  applied  over  a 
horizontal  source  surface  which  moves  toward  a  water  table  have 
been  obtained.  A  commonly  accepted  relationship  between  perme- 
ability and  capillary  pressure  proposed  by  King  (1965)  has  been 
utilized  in  conjunction  with  Darcy's  law  to  formulate  the  mathe- 
matical model.  The  method  of  solution  of  the  model  has  utilized  an 
inverse  formulation  and  finite  difference  techniques.  The  inverse 
formulation  considers  the  cylindrical  coordinates  r  and  z  as  the 
dependent  variables  and  the  potential  function  0  and  the  stream 
function  ijj  as  the  independent  variables  (i.e.  the  problem  is  solved 
for  r  and  z  on  the  $i|j  plane).  Because  of  the  nonlinear  nature  of  the 
equations,  numerical  solutions  are  obtained  by  combining  a 
Newton-Raphson  inner  iteration  with  the  usual  successive  over 
relaxation  outer  iteration.  To  increase  the  accuracy  of  the  numerical 
solution,  a  subdivided  grid  network  has  been  set  up  for  the  region 
between  the  last  two  streamlines  where  the  moisture  content  is 
expected  to  be  extremely  low. 

The  approach  used  for  solving  the  problem  of  partially 
saturated  axisymmetric  infiltration  is  practical  with  modern  digital 
computers.  It  requires  approximately  3  minutes  of  execution  time 
on  a  'UNIVAC  1  108  digital  computer  to  obtain  a  solution  to  a 
problem  such  as  those  presented  herein.  The  computer  output  gives 
the  r  and  z  coordinates  at  each  finite  difference  grid  point  on  the 
$ti|j  plane,  which  can  readily  be  plotted  in  a  flownet  form  to  show 
the  characteristics  of  the  flow  pattern  at  a  glance.  From  the 
computer  solutions,  the  distribution  of  capillary  pressure,  perme- 
ability, or  effective  saturation  over  any  surface  or  plane  of  interest 
can  be  obtained.  Such  quantities  have  been  obtained  from  the 
solutions  and  analyzed  to  define  their  relationship  with  soil 
properties.  The  results  from  these  analyses  are  given  in  a  number  of 
graphs. 

The  solutions  indicate  that  significant  radial  movement  of 
moisture  occurs.  This  "spreading  effect"  is  found  to  be  closely 
related  to  the  hydraulic  properties  of  soils  as  characterized  in  the 
parametric  relationship  used  to  describe  these  properties.  The 
amount  of  radial  movement  increases  as  the  bubbling  pressures  pb 
increases.  It  also  increases  as  the  pore-size  distribution  index  r\ 
increases  and  as  the  other  soil  parameter  b  increases.  The  maximum 
radius  of  the  stream  surface  at  the  water  table,  which  contains  95 
percent  of  the  total  flux,  increases  as  the  value  of  pb  or  b  increases, 
and  decreases  slightly  as  r\  increases. 

While  the  results  are  obtained  based  on  flow  in  porous  media 
concepts,  they  give  insight  into  such  physical  flow  occurrences  in 
nature  as  the  distribution  of  the  permeability  or  effective  saturation 
on  the  surface,  along  the  axis  of  symmetry,  or  on  any  plane 
including  the  axis  of  symmetry  and  how  these  distributions  are 


related  to  the  soil  parameters  defining  the  hydraulic  properties  of 
the  soil.  The  results  show  how  the  moisture  content  decreases  on 
the  surface  beyond  the  area  of  application,  and  how  radial 
movement  into  this  region  causes  a  high  infiltration  rate  at  the  edge 
of  the  source  circle. 

The  following  conclusions  have  been  drawn  directly  from  the 
presentation  and  discussion  of  results  sections.  (The  common 
specifications  of  the  problem  are  R0  =  10',  D  =  16' ,  and  pc/7  = 
1.01  on  the  source  surface.) 

1.  The  "spreading  effect"  due  to  the  radial  movement  of 
moisture  has  been  shown  to  cause  a  higher  infiltration  rate  to  occur 
for  the  axisymmetric  case  than  for  the  one-dimensional  case.  Within 
the  range  of  the  bubbling  pressure  heads  tested,  i.e.  pb/7  varies 
from  1.5  '  to  3.5',  with  n  =  2.50  and  b  =  1.0,  the  radial  movement 
of  moisture  causes  an  excess  of  infiltration  rate  over  the  one- 
dimensional  infiltration  rate  from  19  percent  to  62  percent;  the  19 
percent  increase  being  associated  with  pb/7  =  1-5 '  and  the  60 
percent  with  pb/7  =  3.5' .  When  values  of  b  vary  from  1.0  to  3.0, 
with  t|  =  2.50  and  pb/7  3.50' ,  the  excess  infiltration  rate  increases 
linearly  from  62  percent  to  as  high  as  105  percent.  The  excess 
infiltration  rate  is  also  related  linearly  to  r\  varying  from  62  percent 
to  69  percent  as  n  increases  from  2.50  to  3.50,  with  b  =  1.0  and 
pb/-y  =  3.50' .  The  results  indicate  that  a  higher  percentage  increase 
in  infiltration  rate  can  be  expected  for  soils  with  larger  values  of 
bubbling  pressure,  pore-size  distribution  index  and  b. 

2.  The  radial  movement  of  moisture  has  been  shown  to  cause 
the  streamlines  to  spread  radially  from  the  source  circle  toward  the 
water  table.  The  maximum  radius  of  the  stream  surface  at  the  water 
table  winch  contains  95  percent  of  the  total  flux  expressed  by  r/R0 
increases  from  2.1  to  3.0  as  pb/7  increases  from  1.5  to  3.5,  with  b  = 
1.0  and  n  =  2.50.  It  decreases  slightly  from  2.6  to  2.4  as  n  increases 
from  2.50  to  3.50,  with  b  =  1.0  and  pb/7  =  3.50',  and  increases 
from  2.56  to  2.92  as  b  increases  from  1.0  to  3.0,  with  n  =  2.50  and 
Pb/7  =  3.50'.  Therefore,  the  results  show  that  the  radial  spreading 
of  stream  surfaces  from  the  source  circle  toward  the  water  table  can 
be  significant  if  large  values  of  pb  and  b  are  specified  for  the  soil. 

3.  The  solutions  indicate  that  the  moisture  content  on  the 
ground  surface  beyond  the  source  circle  decreases  as  the  distance 
from  the  source  circle  increases.  For  instance,  the  moisture  content 
decreases  from  approximately  0.8  -  0.99  at  the  edge  of  the  source 
circle  to  about  0.3  -  0.7  at  r/R  0  =  3.0. 

4.  The  solutions  indicate  that  a  definite  edge  effect  causing  a 
high  infiltration  rate  at  the  edge  of  the  source  circle  exists.  The  edge 
effect  decreases  as  the  values  of  pb  or  b  decrease  and  is  not  sensitive 
to  variations  in  n. 
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APPENDICES 


Appendix  A 

One-dimensional  Steady  Downward  Flow 
To  a  Water  Table 


In  an  earlier  study,  Scott  and  Corey  (1961)  used  the 
long-column  technique  of  Childs  and  Collis-George  (1948)  and  the 
controlled  pressure  method  of  Richards  (1931)  to  determine  the 
distribution  of  capillary  pressure  in  columns  of  both  homogeneous 
and  stratified  sands.  An  equation  was  derived  in  terms  of  dimension- 
less  parameters  to  describe  the  distribution  of  capillary  pressure 
during  steady  downward  flow  and  was  referred  as  the  Scott-Corey 
equation.  The  equation  can  be  integrated  to  obtain  an  integral 
expression  for  the  elevation  above  the  water  table,  z. 


Letting  x  =  p0[q0/(  l-bq0)]  1/r|,    Eq.  A-4  can  be  written  as 
1  -bq    «1/t1 


1  -  bq     /  I        q 
o/   \       no 


1 


dx 


0  1    - 


(A-5) 


To  evaluate  I  dx/(l-x  ),  the  the  term  l/U-x^)  is  expressed  by  a 
convergent  series  for  x  <  1  (Arbhabhirama  and  Kridakorn's 
procedure) 


1   -  x 


ti  2ti  3ri 

1    +    X  '    +    X       '    +    X       '    + 


(A-6) 


r/Y 


i 


dP. 


1  -(q/K) 


(A-l) 


in  which  K  is  the  permeability  for  partially  saturated  flow,  pc  is  the 
capillary  pressure,  and  q  is  the  volume  flux  rate.  Arbhabhirama  and 
Kridakorn  (1968)  used  the  dimensionless  form  of  Gardner's 
equation  proposed  by  King  ( 1 965) 


Integrating  the  above  equation  term  by  term  gives 
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dx 


x  + 


0        1  -  x 


2r|+l 


(2T,+  l)(2n  +  2) 


2x 


3r|+l 


(3ti+1)(3t|  +  2) 


+  .   . 
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(A-2) 


+  b 


to  obtain  the  integration  of  Eq.  A-l  by  letting  b  =  1.0.  A  simplified 
expression  describing  the  capillary  pressure  during  steady  downward 
flow  is 


(A-7) 


the  values  of  the  terms  in  the  last  two  parentheses  are  very  small 
when  the  values  of  r\  are  greater  than  2.  Accordingly,  these  terms 
can  be  neglected.  Nearly  all  media  have  values  of  n  greater  than  2. 
Then,  the  integration  can  be  approximated  by 
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dx  x 

=  x  + 


0         1  -x 


n+  1 
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1    - 


T,  +  l 


In      1  - 


q  p 

o    o 
1  -q 


(A-3) 


x  -  — r~r  In    (1  -  x   ) 


T|+  1 


(A-8) 


A  more  general  form  or  a  modified  Arbhabhirama  and 
Kridakorn's  equation  can  be  obtained  by  substituting  Eq.  A-2  into 
Eq.  A-l  without  giving  any  specific  value  of  b. 


Substituting    Eq.   A-8   and   the  expression   for  x   into   Eq.   A-5,   it 
becomes 


-o  ■  r 


dp. 


0  (1-bq^-q^' 


(A-4) 


o  1  -  bq 


i-Lb    i- v- 


r,+  l 


1  -bq 


(A-9) 
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Eq.    A-9   is  a   general   equation   describing   the  capillary   pressure  <i  1 

distribution  during  steady  downward  flow.  K  n, 

°  /P, 


iPh 

For  the  homogeneous  soil,  at  the  region  close  to  the  surface  of 


+  b 


the    soil    column    where    the    capillary    pressure   is   uniform,   and 

•-[((loP./'J/d-t'qo'l   =  0,  or  p0  =  [(l-bq0)/q0]1/rl    or,  in  general  Eq.   A-10   is   used    to   estimate    the   infiltration   rate   of  the  one- 

form  dimensional  steady  downward  infiltration. 
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